1.1. Data manipulation on regional composition (Kenting)

For each site, three types of data will be arranged. The first type is a composition data based on initial categories; the second type is a composition data based on simplified categories; the third type is a data of substrate coverages.

#Get compositions in Kenting (KT)
rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment/Kenting")
library(readxl)
#Import data of Tiaoshih
TSdata<-read_xlsx('TS_composition.xlsx') #Biotic benthic cover
TSdata<-TSdata[,c(3,4,5,7,8)]
TSdata2<-read_xlsx('TS_composition.xlsx', sheet=2) #Substrate cover
TSdata2<-TSdata2[,c(3,4,6,8)]
TSdata<-as.data.frame(TSdata)
TSdata2<-as.data.frame(TSdata2)
TSdata.sum<-aggregate(TSdata[,3],by=list(Mfgroup=TSdata$Mfgroup), FUN=sum )
TSdata.sum.MFG2<-aggregate(TSdata[,3],
                      by=list(Mfgroup=TSdata$MFG2), FUN=sum )
prop<-TSdata.sum[,2]/TSdata2[1,3] #Get coverage
prop2<-TSdata.sum.MFG2[,2]/TSdata2[1,3]#Get coverage
TSdata.sum<-cbind(TSdata.sum,prop)
TSdata.sum.MFG2<-cbind(TSdata.sum.MFG2,prop2)
TSdata.sum[4]<-rep('TS', nrow(TSdata.sum))
TSdata.sum.MFG2[4]<-rep('TS', nrow(TSdata.sum.MFG2))
colnames(TSdata.sum)[4]<-'site'
colnames(TSdata.sum.MFG2)[4]<-'site'
write.csv(TSdata.sum, 'TS_sum.csv') #save table
write.csv(TSdata.sum.MFG2, 'TS_sum_MFG2.csv') #save table
TS2.sum<-aggregate(TSdata2[,3], by=list(Mfgroup=TSdata2$Mfgroup), FUN=sum )#Aggregate substrate coverage
prop<-TS2.sum[,2]/TSdata2[1,3]
TS2.sum$prop<-prop
TS2.sum$site<-rep('TS', nrow(TS2.sum))
write.csv(TS2.sum,'TS_substrates_sum.csv')#Save substrate coverage
###################################################
#Import data of Jialeshui
JLSdata<-read_xlsx('JLS_composition.xlsx')
JLSdata<-JLSdata[,c(2,3,4,6,7)]
JLSdata2<-read_xlsx('JLS_composition.xlsx', sheet=2)
JLSdata2<-JLSdata2[,c(3,4,6,8)]
JLSdata<-as.data.frame(JLSdata)
JLSdata2<-as.data.frame(JLSdata2)
JLSdata.sum<-aggregate(JLSdata[,3],by=list(Mfgroup=JLSdata$Mfgroup), FUN=sum )
JLSdata.sum.MFG2<-aggregate(JLSdata[,3],by=list(Mfgroup=JLSdata$MFG2), FUN=sum )
prop<-JLSdata.sum[,2]/JLSdata2[1,3]
prop2<-JLSdata.sum.MFG2[,2]/JLSdata2[1,3]
JLSdata.sum<-cbind(JLSdata.sum, prop)
JLSdata.sum.MFG2<-cbind(JLSdata.sum.MFG2, prop2)

JLSdata.sum[4]<-rep('JLS', nrow(JLSdata.sum))
JLSdata.sum.MFG2[4]<-rep('JLS', nrow(JLSdata.sum.MFG2))
colnames(JLSdata.sum)[4]<-'site'
colnames(JLSdata.sum.MFG2)[4]<-'site'
write.csv(JLSdata.sum, 'JLS_sum.csv') #save table
write.csv(JLSdata.sum.MFG2, 'JLS_sum_MFG2.csv') #save table
#Aggregate substrate coverage
JLS2.sum<-aggregate(JLSdata2[,3],
                   by=list(Mfgroup=JLSdata2$Mfgroup), FUN=sum )
prop<-JLS2.sum[,2]/JLSdata2[1,3]
JLS2.sum$prop<-prop
JLS2.sum$site<-rep('JLS', nrow(JLS2.sum))
write.csv(JLS2.sum,'JLS_substrates_sum.csv')
###############################################
#Import data of Houbihu
HBHdata<-read_xlsx('HBH_composition.xlsx')
HBHdata<-HBHdata[,c(2,3,4,6,7)]
HBHdata2<-read_xlsx('HBH_composition.xlsx', sheet=2)
HBHdata2<-HBHdata2[,c(3,4,6,8)]
HBHdata<-as.data.frame(HBHdata)
HBHdata2<-as.data.frame(HBHdata2)
HBHdata.sum<-aggregate(HBHdata[,3],by=list(Mfgroup=HBHdata$Mfgroup), FUN=sum )
HBHdata.sum.MFG2<-aggregate(HBHdata[,3],by=list(Mfgroup=HBHdata$MFG2), FUN=sum )
prop<-HBHdata.sum[,2]/HBHdata2[1,3]
prop2<-HBHdata.sum.MFG2[,2]/HBHdata2[1,3]
HBHdata.sum<-cbind(HBHdata.sum, prop)
HBHdata.sum.MFG2<-cbind(HBHdata.sum.MFG2, prop2)
HBHdata.sum<-as.data.frame(HBHdata.sum)
HBHdata.sum[4]<-rep('HBH', nrow(HBHdata.sum))
HBHdata.sum.MFG2[4]<-rep('HBH', nrow(HBHdata.sum.MFG2))
colnames(HBHdata.sum)[4]<-'site'
colnames(HBHdata.sum.MFG2)[4]<-'site'
write.csv(HBHdata.sum, 'HBH_sum.csv') #save table
write.csv(HBHdata.sum.MFG2, 'HBH_sum_MFG2.csv') #save table
#Aggregate substrate coverage
HBH2.sum<-aggregate(HBHdata2[,3],
                    by=list(Mfgroup=HBHdata2$Mfgroup), FUN=sum )
prop<-HBH2.sum[,2]/HBHdata2[1,3]
HBH2.sum$prop<-prop
HBH2.sum$site<-rep('HBH', nrow(HBH2.sum))
write.csv(HBH2.sum,'HBH_substrates_sum.csv')
################################################
#Import data of Leidashih
LDSdata<-read_xlsx('LDS_composition.xlsx')
LDSdata<-LDSdata[,c(3,4,5,7,8)]
LDSdata2<-read_xlsx('LDS_composition.xlsx', sheet=2)
LDSdata2<-LDSdata2[,c(3,4,6)]
LDSdata<-as.data.frame(LDSdata)
LDSdata2<-as.data.frame(LDSdata2)
LDSdata.sum<-aggregate(LDSdata[,3],
        by=list(Mfgroup=LDSdata$Mfgroup), FUN=sum )
LDSdata.sum.MFG2<-aggregate(LDSdata[,3],
                       by=list(Mfgroup=LDSdata$MFG2), FUN=sum )
prop<-LDSdata.sum[,2]/LDSdata2[1,2]
prop2<-LDSdata.sum.MFG2[,2]/LDSdata2[1,2]
LDSdata.sum<-cbind(LDSdata.sum,prop)
LDSdata.sum.MFG2<-cbind(LDSdata.sum.MFG2,prop2)
LDSdata.sum[4]<-rep('LDS', nrow(LDSdata.sum))
LDSdata.sum.MFG2[4]<-rep('LDS', nrow(LDSdata.sum.MFG2))
colnames(LDSdata.sum)[4]<-'site'
colnames(LDSdata.sum.MFG2)[4]<-'site'
write.csv(LDSdata.sum, 'LDS_sum.csv') #save table
write.csv(LDSdata.sum.MFG2, 'LDS_sum_MFG2.csv') #save table
#Aggregate substrate coverage
LDS2.sum<-aggregate(LDSdata2[,2],
                    by=list(Mfgroup=LDSdata2$Mfgroup), FUN=sum )
prop<-LDS2.sum[,2]/LDSdata2[1,2]
LDS2.sum$prop<-prop
LDS2.sum$site<-rep('LDS', nrow(LDS2.sum))
write.csv(LDS2.sum,'LDS_substrates_sum.csv')
####################################################
#Import data of Dinbaisha
DiBSdata<-read_xlsx('DiBS_composition.xlsx')
DiBSdata<-DiBSdata[,c(3,4,5,7,8)]
DiBSdata2<-read_xlsx('DiBS_composition.xlsx', sheet=2)
DiBSdata2<-DiBSdata2[,c(3,4,6,8)]
DiBSdata<-as.data.frame(DiBSdata)
DiBSdata2<-as.data.frame(DiBSdata2)

DiBSdata.sum<-aggregate(DiBSdata[,3],
by=list(Mfgroup=DiBSdata$Mfgroup), FUN=sum )
DiBSdata.sum.MFG2<-aggregate(DiBSdata[,3],
   by=list(Mfgroup=DiBSdata$MFG2), FUN=sum )

prop<-DiBSdata.sum[,2]/DiBSdata2[1,3]
prop2<-DiBSdata.sum.MFG2[,2]/DiBSdata2[1,3]

DiBSdata.sum<-cbind(DiBSdata.sum, prop)
DiBSdata.sum.MFG2<-cbind(DiBSdata.sum.MFG2, prop2)
DiBSdata.sum[4]<-rep('DiBS', nrow(DiBSdata.sum))
DiBSdata.sum.MFG2[4]<-rep('DiBS', nrow(DiBSdata.sum.MFG2))
colnames(DiBSdata.sum)[4]<-'site'
colnames(DiBSdata.sum.MFG2)[4]<-'site'
write.csv(DiBSdata.sum, 'DiBS_sum.csv') #save table
write.csv(DiBSdata.sum.MFG2, 'DiBS_sum_MFG2.csv') #save table

#Aggregate substrate coverage
DiBS2.sum<-aggregate(DiBSdata2[,3],
                    by=list(Mfgroup=DiBSdata2$Mfgroup), FUN=sum )
prop<-DiBS2.sum[,2]/DiBSdata2[1,3]
DiBS2.sum$prop<-prop
DiBS2.sum$site<-rep('DiBS', nrow(DiBS2.sum))
write.csv(DiBS2.sum,'DiBS_substrates_sum.csv')

1.2. Data manipulation on regional composition (Lanyu)

#####Organizing LY data##############
rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment/Lanyu")
library(readxl)
###Import data of GreenGrassland
CCG<-read_xlsx('CCG_composition.xlsx')
CCG<-CCG[,c(3,4,6,8,9)]
CCG2<-read_xlsx('CCG_composition.xlsx', sheet=2)
CCG2<-CCG2[,c(3,4,6,8)]
CCG<-as.data.frame(CCG)
CCG2<-as.data.frame(CCG2)

CCG.sum<-aggregate(CCG[,3],
                        by=list(Mfgroup=CCG$Mfgroup), FUN=sum )
CCG.MFG2.sum<-aggregate(CCG[,3],
                  by=list(Mfgroup=CCG$MFG2), FUN=sum )

prop<-CCG.sum[,2]/CCG2[1,3]
prop<-as.data.frame(prop)

CCG.sum<-cbind(CCG.sum, prop)
CCG.sum[4]<-rep('CCG', nrow(CCG.sum))
colnames(CCG.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(CCG.sum,'CCG_sum.csv')

prop<-CCG.MFG2.sum[,2]/CCG2[1,3]
prop<-as.data.frame(prop)

CCG.MFG2.sum<-cbind(CCG.MFG2.sum, prop)
CCG.MFG2.sum[4]<-rep('CCG', nrow(CCG.MFG2.sum))
colnames(CCG.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(CCG.MFG2.sum,'CCG_MFG2_sum.csv')

#Substrate coverages
CCG2.sum<-aggregate(CCG2[,3],
          by=list(Mfgroup=CCG2$Mfgroup), FUN=sum )
prop<-CCG2.sum[,2]/CCG2[1,3]
CCG2.sum$prop<-prop
CCG2.sum$site<-rep('CCG', nrow(CCG2.sum))
write.csv(CCG2.sum,'CCG_substrates_sum.csv')

######Import data of HenRock#
HR<-read_xlsx('HR_composition.xlsx')
HR<-HR[,c(3,4,6,8,9)]
HR2<-read_xlsx('HR_composition.xlsx', sheet=2)
HR2<-HR2[,c(3,4,6,8)]
HR<-as.data.frame(HR)
HR2<-as.data.frame(HR2)

HR.sum<-aggregate(HR[,3],
                   by=list(Mfgroup=HR$Mfgroup), FUN=sum )
HR.MFG2.sum<-aggregate(HR[,3],
                        by=list(Mfgroup=HR$MFG2), FUN=sum )

prop<-HR.sum[,2]/HR2[1,3]
prop<-as.data.frame(prop)

HR.sum<-cbind(HR.sum, prop)
HR.sum[4]<-rep('HR', nrow(HR.sum))
colnames(HR.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(HR.sum,'HR_sum.csv')

prop<-HR.MFG2.sum[,2]/HR2[1,3]
prop<-as.data.frame(prop)

HR.MFG2.sum<-cbind(HR.MFG2.sum, prop)
HR.MFG2.sum[4]<-rep('HR', nrow(HR.MFG2.sum))
colnames(HR.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(HR.MFG2.sum,'HR_MFG2_sum.csv')

#substrates coverage
HR2.sum<-aggregate(HR2[,3],
                    by=list(Mfgroup=HR2$Mfgroup), FUN=sum )
prop<-HR2.sum[,2]/HR2[1,3]
HR2.sum$prop<-prop
HR2.sum$site<-rep('HR', nrow(HR2.sum))
write.csv(HR2.sum,'HR_substrates_sum.csv')

##Import data of Dongchin#
DC<-read_xlsx('DC_composition.xlsx')
DC<-DC[,c(3,4,6,8,9)]
DC2<-read_xlsx('DC_composition.xlsx', sheet=2)
DC2<-DC2[,c(3,4,6,8)]
DC<-as.data.frame(DC)
DC2<-as.data.frame(DC2)

DC.sum<-aggregate(DC[,3],
                  by=list(Mfgroup=DC$Mfgroup), FUN=sum )
DC.MFG2.sum<-aggregate(DC[,3],
                       by=list(Mfgroup=DC$MFG2), FUN=sum )

prop<-DC.sum[,2]/DC2[1,3]
prop<-as.data.frame(prop)

DC.sum<-cbind(DC.sum, prop)
DC.sum[4]<-rep('DC', nrow(DC.sum))
colnames(DC.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(DC.sum,'DC_sum.csv')

prop<-DC.MFG2.sum[,2]/DC2[1,3]
prop<-as.data.frame(prop)

DC.MFG2.sum<-cbind(DC.MFG2.sum, prop)
DC.MFG2.sum[4]<-rep('DC', nrow(DC.MFG2.sum))
colnames(DC.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(DC.MFG2.sum,'DC_MFG2_sum.csv')

#substrate coverage#
DC2.sum<-aggregate(DC2[,3],
                   by=list(Mfgroup=DC2$Mfgroup), FUN=sum )
prop<-DC2.sum[,2]/DC2[1,3]
DC2.sum$prop<-prop
DC2.sum$site<-rep('DC', nrow(DC2.sum))
write.csv(DC2.sum,'DC_substrates_sum.csv')

##Import data of TwoLions#
TL<-read_xlsx('TL_composition.xlsx')
TL<-TL[,c(3,4,6,8,9)]
TL2<-read_xlsx('TL_composition.xlsx', sheet=2)
TL2<-TL2[,c(3,4,6,8)]
TL<-as.data.frame(TL)
TL2<-as.data.frame(TL2)

TL.sum<-aggregate(TL[,3],
                  by=list(Mfgroup=TL$Mfgroup), FUN=sum )
TL.MFG2.sum<-aggregate(TL[,3],
                       by=list(Mfgroup=TL$MFG2), FUN=sum )

prop<-TL.sum[,2]/TL2[1,3]
prop<-as.data.frame(prop)

TL.sum<-cbind(TL.sum, prop)
TL.sum[4]<-rep('TL', nrow(TL.sum))
colnames(TL.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(TL.sum,'TL_sum.csv')

prop<-TL.MFG2.sum[,2]/TL2[1,3]
prop<-as.data.frame(prop)

TL.MFG2.sum<-cbind(TL.MFG2.sum, prop)
TL.MFG2.sum[4]<-rep('TL', nrow(TL.MFG2.sum))
colnames(TL.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(TL.MFG2.sum,'TL_MFG2_sum.csv')

#substrates coverage#
TL2.sum<-aggregate(TL2[,3],
                   by=list(Mfgroup=TL2$Mfgroup), FUN=sum )
prop<-TL2.sum[,2]/TL2[1,3]
TL2.sum$prop<-prop
TL2.sum$site<-rep('TL', nrow(TL2.sum))
write.csv(TL2.sum,'TL_substrates_sum.csv')

##Import data of YaYo#
YY<-read_xlsx('YY_composition.xlsx')
YY<-YY[,c(3,4,6,8,9)]
YY2<-read_xlsx('YY_composition.xlsx', sheet=2)
YY2<-YY2[,c(3,4,6,8)]
YY<-as.data.frame(YY)
YY2<-as.data.frame(YY2)

YY.sum<-aggregate(YY[,3],
                  by=list(Mfgroup=YY$Mfgroup), FUN=sum )
YY.MFG2.sum<-aggregate(YY[,3],
                       by=list(Mfgroup=YY$MFG2), FUN=sum )

prop<-YY.sum[,2]/YY2[1,3]
prop<-as.data.frame(prop)

YY.sum<-cbind(YY.sum, prop)
YY.sum[4]<-rep('YY', nrow(YY.sum))
colnames(YY.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(YY.sum,'YY_sum.csv')


prop<-YY.MFG2.sum[,2]/YY2[1,3]
prop<-as.data.frame(prop)

YY.MFG2.sum<-cbind(YY.MFG2.sum, prop)
YY.MFG2.sum[4]<-rep('YY', nrow(YY.MFG2.sum))
colnames(YY.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(YY.MFG2.sum,'YY_MFG2_sum.csv')

#substrate coverage
YY2.sum<-aggregate(YY2[,3],
                   by=list(Mfgroup=YY2$Mfgroup), FUN=sum )
prop<-YY2.sum[,2]/YY2[1,3]
YY2.sum$prop<-prop
YY2.sum$site<-rep('YY', nrow(YY2.sum))
write.csv(YY2.sum,'YY_substrates_sum.csv')

1.3. Data manipulation on regional composition (Ludao)

#####Organizing Ludao data##############
rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment/Ludao")
library(readxl)

##Import data of Guiwan
GW<-read_xlsx('GW_composition.xlsx')
GW<-GW[,c(3,4,6,8,9)]
GW2<-read_xlsx('GW_composition.xlsx', sheet=2)
GW2<-GW2[,c(3,4,6,8)]
GW<-as.data.frame(GW)
GW2<-as.data.frame(GW2)

GW.sum<-aggregate(GW[,3],
                        by=list(Mfgroup=GW$Mfgroup), FUN=sum )
GW.MFG2.sum<-aggregate(GW[,3],
                  by=list(Mfgroup=GW$MFG2), FUN=sum )

prop<-GW.sum[,2]/GW2[1,3]
prop<-as.data.frame(prop)

GW.sum<-cbind(GW.sum, prop)
GW.sum[4]<-rep('GW', nrow(GW.sum))
colnames(GW.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(GW.sum,'GW_sum.csv')

prop<-GW.MFG2.sum[,2]/GW2[1,3]
prop<-as.data.frame(prop)

GW.MFG2.sum<-cbind(GW.MFG2.sum, prop)
GW.MFG2.sum[4]<-rep('GW', nrow(GW.MFG2.sum))
colnames(GW.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(GW.MFG2.sum,'GW_MFG2_sum.csv')

#Substrate coverage
GW2.sum<-aggregate(GW2[,3],
          by=list(Mfgroup=GW2$Mfgroup), FUN=sum )
prop<-GW2.sum[,2]/GW2[1,3]
GW2.sum$prop<-prop
GW2.sum$site<-rep('GW', nrow(GW2.sum))
write.csv(GW2.sum,'GW_substrates_sum.csv')


##Import data of Chaikou
CK<-read_xlsx('CK_composition.xlsx')
CK<-CK[,c(3,4,6,8,9)]
CK2<-read_xlsx('CK_composition.xlsx', sheet=2)
CK2<-CK2[,c(3,4,6,8)]
CK<-as.data.frame(CK)
CK2<-as.data.frame(CK2)

CK.sum<-aggregate(CK[,3],
                  by=list(Mfgroup=CK$Mfgroup), FUN=sum )
CK.MFG2.sum<-aggregate(CK[,3],
                       by=list(Mfgroup=CK$MFG2), FUN=sum )

prop<-CK.sum[,2]/CK2[1,3]
prop<-as.data.frame(prop)

CK.sum<-cbind(CK.sum, prop)
CK.sum[4]<-rep('CK', nrow(CK.sum))
colnames(CK.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(CK.sum,'CK_sum.csv')

prop<-CK.MFG2.sum[,2]/CK2[1,3]
prop<-as.data.frame(prop)

CK.MFG2.sum<-cbind(CK.MFG2.sum, prop)
CK.MFG2.sum[4]<-rep('CK', nrow(CK.MFG2.sum))
colnames(CK.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(CK.MFG2.sum,'CK_MFG2_sum.csv')


#substrate coverage
CK2.sum<-aggregate(CK2[,3],
                   by=list(Mfgroup=CK2$Mfgroup), FUN=sum )
prop<-CK2.sum[,2]/CK2[1,3]
CK2.sum$prop<-prop
CK2.sum$site<-rep('CK', nrow(CK2.sum))
write.csv(CK2.sum,'CK_substrates_sum.csv')


##Import data of Shihlang#
SL<-read_xlsx('SL_composition.xlsx')
SL<-SL[,c(3,4,6,8,9)]
SL2<-read_xlsx('SL_composition.xlsx', sheet=2)
SL2<-SL2[,c(3,4,6,8)]
SL<-as.data.frame(SL)
SL2<-as.data.frame(SL2)

SL.sum<-aggregate(SL[,3],
                  by=list(Mfgroup=SL$Mfgroup), FUN=sum )
SL.MFG2.sum<-aggregate(SL[,3],
                       by=list(Mfgroup=SL$MFG2), FUN=sum )

prop<-SL.sum[,2]/SL2[1,3]
prop<-as.data.frame(prop)

SL.sum<-cbind(SL.sum, prop)
SL.sum[4]<-rep('SL', nrow(SL.sum))
colnames(SL.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(SL.sum,'SL_sum.csv')

prop<-SL.MFG2.sum[,2]/SL2[1,3]
prop<-as.data.frame(prop)

SL.MFG2.sum<-cbind(SL.MFG2.sum, prop)
SL.MFG2.sum[4]<-rep('SL', nrow(SL.MFG2.sum))
colnames(SL.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(SL.MFG2.sum,'SL_MFG2_sum.csv')


#Substrate coverage
SL2.sum<-aggregate(SL2[,3],
                   by=list(Mfgroup=SL2$Mfgroup), FUN=sum )
prop<-SL2.sum[,2]/SL2[1,3]
SL2.sum$prop<-prop
SL2.sum$site<-rep('SL', nrow(SL2.sum))
write.csv(SL2.sum,'SL_substrates_sum.csv')

##Import Dabaisha data
DaBS<-read_xlsx('DaBS_composition.xlsx')
DaBS<-DaBS[,c(3,4,6,8,9)]
DaBS2<-read_xlsx('DaBS_composition.xlsx', sheet=2)
DaBS2<-DaBS2[,c(3,4,6,8)]
DaBS<-as.data.frame(DaBS)
DaBS2<-as.data.frame(DaBS2)

DaBS.sum<-aggregate(DaBS[,3],
                  by=list(Mfgroup=DaBS$Mfgroup), FUN=sum )
DaBS.MFG2.sum<-aggregate(DaBS[,3],
                       by=list(Mfgroup=DaBS$MFG2), FUN=sum )

prop<-DaBS.sum[,2]/DaBS2[1,3]
prop<-as.data.frame(prop)

DaBS.sum<-cbind(DaBS.sum, prop)
DaBS.sum[4]<-rep('DaBS', nrow(DaBS.sum))
colnames(DaBS.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(DaBS.sum,'DaBS_sum.csv')


prop<-DaBS.MFG2.sum[,2]/DaBS2[1,3]
prop<-as.data.frame(prop)

DaBS.MFG2.sum<-cbind(DaBS.MFG2.sum, prop)
DaBS.MFG2.sum[4]<-rep('DaBS', nrow(DaBS.MFG2.sum))
colnames(DaBS.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(DaBS.MFG2.sum,'DaBS_MFG2_sum.csv')

#Substrate coverage
DaBS2.sum<-aggregate(DaBS2[,3],
                   by=list(Mfgroup=DaBS2$Mfgroup), FUN=sum )
prop<-DaBS2.sum[,2]/DaBS2[1,3]
DaBS2.sum$prop<-prop
DaBS2.sum$site<-rep('DaBS', nrow(DaBS2.sum))
write.csv(DaBS2.sum,'DaBS_substrates_sum.csv')

##Import data of Gongguanbi#
GGB<-read_xlsx('GGB_composition.xlsx')
GGB<-GGB[,c(3,4,6,8,9)]
GGB2<-read_xlsx('GGB_composition.xlsx', sheet=2)
GGB2<-GGB2[,c(3,4,6,8)]
GGB<-as.data.frame(GGB)
GGB2<-as.data.frame(GGB2)

GGB.sum<-aggregate(GGB[,3],
                  by=list(Mfgroup=GGB$Mfgroup), FUN=sum )
GGB.MFG2.sum<-aggregate(GGB[,3],
                       by=list(Mfgroup=GGB$MFG2), FUN=sum )

prop<-GGB.sum[,2]/GGB2[1,3]
prop<-as.data.frame(prop)

GGB.sum<-cbind(GGB.sum, prop)
GGB.sum[4]<-rep('GGB', nrow(GGB.sum))
colnames(GGB.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(GGB.sum,'GGB_sum.csv')


prop<-GGB.MFG2.sum[,2]/GGB2[1,3]
prop<-as.data.frame(prop)

GGB.MFG2.sum<-cbind(GGB.MFG2.sum, prop)
GGB.MFG2.sum[4]<-rep('GGB', nrow(GGB.MFG2.sum))
colnames(GGB.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(GGB.MFG2.sum,'GGB_MFG2_sum.csv')

#Substrate coverage
GGB2.sum<-aggregate(GGB2[,3],
                   by=list(Mfgroup=GGB2$Mfgroup), FUN=sum )
prop<-GGB2.sum[,2]/GGB2[1,3]
GGB2.sum$prop<-prop
GGB2.sum$site<-rep('GGB', nrow(GGB2.sum))
write.csv(GGB2.sum,'GGB_substrates_sum.csv')

1.4. Data manipulation on regional composition (East Coast)

rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment/East_Coast")
library(readxl)
##Import data of Jihui
JH<-read_xlsx('JH_composition.xlsx')
JH<-JH[,c(3,4,6,8,9)]
JH2<-read_xlsx('JH_composition.xlsx', sheet=2)
JH2<-JH2[,c(3,4,6,8)]
JH<-as.data.frame(JH)
JH2<-as.data.frame(JH2)

JH.sum<-aggregate(JH[,3],
                   by=list(Mfgroup=JH$Mfgroup), FUN=sum )
JH.MFG2.sum<-aggregate(JH[,3],
                        by=list(Mfgroup=JH$MFG2), FUN=sum )

prop<-JH.sum[,2]/JH2[1,3]
prop<-as.data.frame(prop)

JH.sum<-cbind(JH.sum, prop)
JH.sum[4]<-rep('JH', nrow(JH.sum))
colnames(JH.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(JH.sum,'JH_sum.csv')


prop<-JH.MFG2.sum[,2]/JH2[1,3]
prop<-as.data.frame(prop)

JH.MFG2.sum<-cbind(JH.MFG2.sum, prop)
JH.MFG2.sum[4]<-rep('JH', nrow(JH.MFG2.sum))
colnames(JH.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(JH.MFG2.sum,'JH_MFG2_sum.csv')

#Substrate coverage
JH2.sum<-aggregate(JH2[,3],
                    by=list(Mfgroup=JH2$Mfgroup), FUN=sum )
prop<-JH2.sum[,2]/JH2[1,3]
JH2.sum$prop<-prop
JH2.sum$site<-rep('JH', nrow(JH2.sum))
write.csv(JH2.sum,'JH_substrates_sum.csv')


##Import data of Jiamuziwan#
JMZ<-read_xlsx('JMZ_composition.xlsx')
JMZ<-JMZ[,c(3,4,6,8,9)]
JMZ2<-read_xlsx('JMZ_composition.xlsx', sheet=2)
JMZ2<-JMZ2[,c(3,4,6,8)]
JMZ<-as.data.frame(JMZ)
JMZ2<-as.data.frame(JMZ2)

JMZ.sum<-aggregate(JMZ[,3],
                  by=list(Mfgroup=JMZ$Mfgroup), FUN=sum )
JMZ.MFG2.sum<-aggregate(JMZ[,3],
                       by=list(Mfgroup=JMZ$MFG2), FUN=sum )

prop<-JMZ.sum[,2]/JMZ2[1,3]
prop<-as.data.frame(prop)

JMZ.sum<-cbind(JMZ.sum, prop)
JMZ.sum[4]<-rep('JMZ', nrow(JMZ.sum))
colnames(JMZ.sum)[c(2,4)]<-c('area','site')

#export table based on initial categories
write.csv(JMZ.sum,'JMZ_sum.csv')

prop<-JMZ.MFG2.sum[,2]/JMZ2[1,3]
prop<-as.data.frame(prop)
JMZ.MFG2.sum<-cbind(JMZ.MFG2.sum, prop)
JMZ.MFG2.sum[4]<-rep('JMZ', nrow(JMZ.MFG2.sum))
colnames(JMZ.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(JMZ.MFG2.sum,'JMZ_MFG2_sum.csv')


#substrate coverage
JMZ2.sum<-aggregate(JMZ2[,3],
                   by=list(Mfgroup=JMZ2$Mfgroup), FUN=sum )
prop<-JMZ2.sum[,2]/JMZ2[1,3]
JMZ2.sum$prop<-prop
JMZ2.sum$site<-rep('JMZ', nrow(JMZ2.sum))
write.csv(JMZ2.sum,'JMZ_substrates_sum.csv')

##Import data of Fenniaolin#
FNL<-read_xlsx('FNL_composition.xlsx')
FNL<-FNL[,c(3,4,6,8,9)]
FNL2<-read_xlsx('FNL_composition.xlsx', sheet=2)
FNL2<-FNL2[,c(3,4,6,8)]
FNL<-as.data.frame(FNL)
FNL2<-as.data.frame(FNL2)

FNL.sum<-aggregate(FNL[,3],
                   by=list(Mfgroup=FNL$Mfgroup), FUN=sum )
FNL.MFG2.sum<-aggregate(FNL[,3],
                        by=list(Mfgroup=FNL$MFG2), FUN=sum )

prop<-FNL.sum[,2]/FNL2[1,3]
prop<-as.data.frame(prop)

FNL.sum<-cbind(FNL.sum, prop)
FNL.sum[4]<-rep('FNL', nrow(FNL.sum))
colnames(FNL.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(FNL.sum,'FNL_sum.csv')


prop<-FNL.MFG2.sum[,2]/FNL2[1,3]
prop<-as.data.frame(prop)

FNL.MFG2.sum<-cbind(FNL.MFG2.sum, prop)
FNL.MFG2.sum[4]<-rep('FNL', nrow(FNL.MFG2.sum))
colnames(FNL.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(FNL.MFG2.sum,'FNL_MFG2_sum.csv')


#substrate coverage
FNL2.sum<-aggregate(FNL2[,3],
                    by=list(Mfgroup=FNL2$Mfgroup), FUN=sum )
prop<-FNL2.sum[,2]/FNL2[1,3]
FNL2.sum$prop<-prop
FNL2.sum$site<-rep('FNL', nrow(FNL2.sum))
write.csv(FNL2.sum,'FNL_substrates_sum.csv')

##Import data of Hsinshe#
HS<-read_xlsx('HS_composition.xlsx')
HS<-HS[,c(3,4,6,8,9)]
HS2<-read_xlsx('HS_composition.xlsx', sheet=2)
HS2<-HS2[,c(3,4,6,8)]
HS<-as.data.frame(HS)
HS2<-as.data.frame(HS2)

HS.sum<-aggregate(HS[,3],
                   by=list(Mfgroup=HS$Mfgroup), FUN=sum )
HS.MFG2.sum<-aggregate(HS[,3],
                        by=list(Mfgroup=HS$MFG2), FUN=sum )

prop<-HS.sum[,2]/HS2[1,3]
prop<-as.data.frame(prop)

HS.sum<-cbind(HS.sum, prop)
HS.sum[4]<-rep('HS', nrow(HS.sum))
colnames(HS.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(HS.sum,'HS_sum.csv')


prop<-HS.MFG2.sum[,2]/HS2[1,3]
prop<-as.data.frame(prop)

HS.MFG2.sum<-cbind(HS.MFG2.sum, prop)
HS.MFG2.sum[4]<-rep('HS', nrow(HS.MFG2.sum))
colnames(HS.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(HS.MFG2.sum,'HS_MFG2_sum.csv')

#Substrate coverage
HS2.sum<-aggregate(HS2[,3],
                    by=list(Mfgroup=HS2$Mfgroup), FUN=sum )
prop<-HS2.sum[,2]/HS2[1,3]
HS2.sum$prop<-prop
HS2.sum$site<-rep('HS', nrow(HS2.sum))
write.csv(HS2.sum,'HS_substrates_sum.csv')

##Import data of Jiqi
JQ<-read_xlsx('JQ_composition.xlsx')
JQ<-JQ[,c(3,4,6,8,9)]
JQ2<-read_xlsx('JQ_composition.xlsx', sheet=2)
JQ2<-JQ2[,c(3,4,6,8)]
JQ<-as.data.frame(JQ)
JQ2<-as.data.frame(JQ2)

JQ.sum<-aggregate(JQ[,3],
                  by=list(Mfgroup=JQ$Mfgroup), FUN=sum )
JQ.MFG2.sum<-aggregate(JQ[,3],
                       by=list(Mfgroup=JQ$MFG2), FUN=sum )

prop<-JQ.sum[,2]/JQ2[1,3]
prop<-as.data.frame(prop)

JQ.sum<-cbind(JQ.sum, prop)
JQ.sum[4]<-rep('JQ', nrow(JQ.sum))
colnames(JQ.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(JQ.sum,'JQ_sum.csv')


prop<-JQ.MFG2.sum[,2]/JQ2[1,3]
prop<-as.data.frame(prop)

JQ.MFG2.sum<-cbind(JQ.MFG2.sum, prop)
JQ.MFG2.sum[4]<-rep('JQ', nrow(JQ.MFG2.sum))
colnames(JQ.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(JQ.MFG2.sum,'JQ_MFG2_sum.csv')

#substrate coverage
JQ2.sum<-aggregate(JQ2[,3],
                   by=list(Mfgroup=JQ2$Mfgroup), FUN=sum )
prop<-JQ2.sum[,2]/JQ2[1,3]
JQ2.sum$prop<-prop
JQ2.sum$site<-rep('JQ', nrow(JQ2.sum))
write.csv(JQ2.sum,'JQ_substrates_sum.csv')

1.5. Data manipulation on regional composition (North Coast)

rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment/North_Coast")
library(readxl)
##Import data of Bitou#
BT<-read_xlsx('BT_composition.xlsx')
BT<-BT[,c(3,4,6,8,9)]
BT2<-read_xlsx('BT_composition.xlsx', sheet=2)
BT2<-BT2[,c(3,4,6,8)]
BT<-as.data.frame(BT)
BT2<-as.data.frame(BT2)

BT.sum<-aggregate(BT[,3],
                  by=list(Mfgroup=BT$Mfgroup), FUN=sum )
BT.MFG2.sum<-aggregate(BT[,3],
                       by=list(Mfgroup=BT$MFG2), FUN=sum )

prop<-BT.sum[,2]/BT2[1,3]
prop<-as.data.frame(prop)

BT.sum<-cbind(BT.sum, prop)
BT.sum[4]<-rep('BT', nrow(BT.sum))
colnames(BT.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(BT.sum,'BT_sum.csv')


prop<-BT.MFG2.sum[,2]/BT2[1,3]
prop<-as.data.frame(prop)

BT.MFG2.sum<-cbind(BT.MFG2.sum, prop)
BT.MFG2.sum[4]<-rep('BT', nrow(BT.MFG2.sum))
colnames(BT.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(BT.MFG2.sum,'BT_MFG2_sum.csv')

#substrate coverage
BT2.sum<-aggregate(BT2[,3],
                   by=list(Mfgroup=BT2$Mfgroup), FUN=sum )
prop<-BT2.sum[,2]/BT2[1,3]
BT2.sum$prop<-prop
BT2.sum$site<-rep('BT', nrow(BT2.sum))
write.csv(BT2.sum,'BT_substrates_sum.csv')

##Import data of Chaojing#
CJ<-read_xlsx('CJ_composition.xlsx')
CJ<-CJ[,c(3,4,6,8,9)]
CJ2<-read_xlsx('CJ_composition.xlsx', sheet=2)
CJ2<-CJ2[,c(3,4,6,8)]
CJ<-as.data.frame(CJ)
CJ2<-as.data.frame(CJ2)

CJ.sum<-aggregate(CJ[,3],
                  by=list(Mfgroup=CJ$Mfgroup), FUN=sum )
CJ.MFG2.sum<-aggregate(CJ[,3],
                       by=list(Mfgroup=CJ$MFG2), FUN=sum )

prop<-CJ.sum[,2]/CJ2[1,3]
prop<-as.data.frame(prop)

CJ.sum<-cbind(CJ.sum, prop)
CJ.sum[4]<-rep('CJ', nrow(CJ.sum))
colnames(CJ.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(CJ.sum,'CJ_sum.csv')


prop<-CJ.MFG2.sum[,2]/CJ2[1,3]
prop<-as.data.frame(prop)

CJ.MFG2.sum<-cbind(CJ.MFG2.sum, prop)
CJ.MFG2.sum[4]<-rep('CJ', nrow(CJ.MFG2.sum))
colnames(CJ.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(CJ.MFG2.sum,'CJ_MFG2_sum.csv')

#substrate coverage
CJ2.sum<-aggregate(CJ2[,3],
                   by=list(Mfgroup=CJ2$Mfgroup), FUN=sum )
prop<-CJ2.sum[,2]/CJ2[1,3]
CJ2.sum$prop<-prop
CJ2.sum$site<-rep('CJ', nrow(CJ2.sum))
write.csv(CJ2.sum,'CJ_substrates_sum.csv')


##Import data of Longdong
LD<-read_xlsx('LD_composition.xlsx')
LD<-LD[,c(3,4,6,8,9)]
LD2<-read_xlsx('LD_composition.xlsx', sheet=2)
LD2<-LD2[,c(3,4,6,8)]
LD<-as.data.frame(LD)
LD2<-as.data.frame(LD2)

LD.sum<-aggregate(LD[,3],
                  by=list(Mfgroup=LD$Mfgroup), FUN=sum )
LD.MFG2.sum<-aggregate(LD[,3],
                       by=list(Mfgroup=LD$MFG2), FUN=sum )

prop<-LD.sum[,2]/LD2[1,3]
prop<-as.data.frame(prop)

LD.sum<-cbind(LD.sum, prop)
LD.sum[4]<-rep('LD', nrow(LD.sum))
colnames(LD.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(LD.sum,'LD_sum.csv')


prop<-LD.MFG2.sum[,2]/LD2[1,3]
prop<-as.data.frame(prop)

LD.MFG2.sum<-cbind(LD.MFG2.sum, prop)
LD.MFG2.sum[4]<-rep('LD', nrow(LD.MFG2.sum))
colnames(LD.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(LD.MFG2.sum,'LD_MFG2_sum.csv')

#substrate coverage
LD2.sum<-aggregate(LD2[,3],
                   by=list(Mfgroup=LD2$Mfgroup), FUN=sum )
prop<-LD2.sum[,2]/LD2[1,3]
LD2.sum$prop<-prop
LD2.sum$site<-rep('LD', nrow(LD2.sum))
write.csv(LD2.sum,'LD_substrates_sum.csv')

##Import data of 82.5K#
ZM<-read_xlsx('825_composition.xlsx')
ZM<-ZM[,c(3,4,6,8,9)]
ZM2<-read_xlsx('825_composition.xlsx', sheet=2)
ZM2<-ZM2[,c(3,4,6,8)]
ZM<-as.data.frame(ZM)
ZM2<-as.data.frame(ZM2)

ZM.sum<-aggregate(ZM[,3],
                  by=list(Mfgroup=ZM$Mfgroup), FUN=sum )
ZM.MFG2.sum<-aggregate(ZM[,3],
                       by=list(Mfgroup=ZM$MFG2), FUN=sum )

prop<-ZM.sum[,2]/ZM2[1,3]
prop<-as.data.frame(prop)

ZM.sum<-cbind(ZM.sum, prop)
ZM.sum[4]<-rep('82.5', nrow(ZM.sum))
colnames(ZM.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(ZM.sum,'82.5_sum.csv')

prop<-ZM.MFG2.sum[,2]/ZM2[1,3]
prop<-as.data.frame(prop)

ZM.MFG2.sum<-cbind(ZM.MFG2.sum, prop)
ZM.MFG2.sum[4]<-rep('82.5', nrow(ZM.MFG2.sum))
colnames(ZM.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(ZM.MFG2.sum,'82.5_MFG2_sum.csv')


#substrate coverage
ZM2.sum<-aggregate(ZM2[,3],
                   by=list(Mfgroup=ZM2$Mfgroup), FUN=sum )
prop<-ZM2.sum[,2]/ZM2[1,3]
ZM2.sum$prop<-prop
ZM2.sum$site<-rep('82.5', nrow(ZM2.sum))
write.csv(ZM2.sum,'82.5_substrates_sum.csv')


##Import data of Meiyenshan#
MYS<-read_xlsx('MYS_composition.xlsx')
MYS<-MYS[,c(3,4,6,8,9)]
MYS2<-read_xlsx('MYS_composition.xlsx', sheet=2)
MYS2<-MYS2[,c(3,4,6,8)]
MYS<-as.data.frame(MYS)
MYS2<-as.data.frame(MYS2)

MYS.sum<-aggregate(MYS[,3],
                  by=list(Mfgroup=MYS$Mfgroup), FUN=sum )
MYS.MFG2.sum<-aggregate(MYS[,3],
                       by=list(Mfgroup=MYS$MFG2), FUN=sum )

prop<-MYS.sum[,2]/MYS2[1,3]
prop<-as.data.frame(prop)

MYS.sum<-cbind(MYS.sum, prop)
MYS.sum[4]<-rep('MYS', nrow(MYS.sum))
colnames(MYS.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(MYS.sum,'MYS_sum.csv')


prop<-MYS.MFG2.sum[,2]/MYS2[1,3]
prop<-as.data.frame(prop)

MYS.MFG2.sum<-cbind(MYS.MFG2.sum, prop)
MYS.MFG2.sum[4]<-rep('MYS', nrow(MYS.MFG2.sum))
colnames(MYS.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(MYS.MFG2.sum,'MYS_MFG2_sum.csv')

#substrate coverage
MYS2.sum<-aggregate(MYS2[,3],
                   by=list(Mfgroup=MYS2$Mfgroup), FUN=sum )
prop<-MYS2.sum[,2]/MYS2[1,3]
MYS2.sum$prop<-prop
MYS2.sum$site<-rep('MYS', nrow(BT2.sum))
write.csv(MYS2.sum,'MYS_substrates_sum.csv')

2. Data manipulation to combine all regional compositions.

This step is to combine all composition data with simplified categories, forming a big data frame. Within this data frame, each row tells one particular site, and each column tells one particular category.

rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment")
######Import all data based on simplified categories
TS<-read.csv('~/myExperiment/Kenting/TS_sum_MFG2.csv')
LDS<-read.csv('~/myExperiment/Kenting/LDS_sum_MFG2.csv')
JLS<-read.csv('~/myExperiment/Kenting/JLS_sum_MFG2.csv')
HBH<-read.csv('~/myExperiment/Kenting/HBH_sum_MFG2.csv')
DiBS<-read.csv('~/myExperiment/Kenting/DiBS_sum_MFG2.csv')

GW<-read.csv('~/myExperiment/Ludao/GW_MFG2_sum.csv')
CK<-read.csv('~/myExperiment/Ludao/CK_MFG2_sum.csv')
SL<-read.csv('~/myExperiment/Ludao/SL_MFG2_sum.csv')
DaBS<-read.csv('~/myExperiment/Ludao/DaBS_MFG2_sum.csv')
GGB<-read.csv('~/myExperiment/Ludao/GGB_MFG2_sum.csv')

CCG<-read.csv('~/myExperiment/Lanyu/CCG_MFG2_sum.csv')
YY<-read.csv('~/myExperiment/Lanyu/YY_MFG2_sum.csv')
HR<-read.csv('~/myExperiment/Lanyu/HR_MFG2_sum.csv')
TL<-read.csv('~/myExperiment/Lanyu/TL_MFG2_sum.csv')
DC<-read.csv('~/myExperiment/Lanyu/DC_MFG2_Sum.csv')

JH<-read.csv('~/myExperiment/East Coast/JH_MFG2_sum.csv')
JMZ<-read.csv('~/myExperiment/East Coast/JMZ_MFG2_sum.csv')
HS<-read.csv('~/myExperiment/East Coast/HS_MFG2_sum.csv')
JQ<-read.csv('~/myExperiment/East Coast/JQ_MFG2_sum.csv')
FNL<-read.csv('~/myExperiment/East Coast/FNL_MFG2_sum.csv')

LD<-read.csv('~/myExperiment/North Coast/LD_MFG2_sum.csv')
BT<-read.csv('~/myExperiment/North Coast/BT_MFG2_sum.csv')
CJ<-read.csv('~/myExperiment/North Coast/CJ_MFG2_sum.csv')
MYS<-read.csv('~/myExperiment/North Coast/MYS_MFG2_sum.csv')
ZM<-read.csv('~/myExperiment/North Coast/82.5_MFG2_sum.csv')

###Organizing the datasets###
#Set the factor telling regions
colnames(TS)[3]<-c('area')
colnames(DiBS)[3]<-c('area')
colnames(JLS)[3]<-c('area')
colnames(LDS)[3]<-c('area')
colnames(HBH)[3]<-c('area')

colnames(TS)[4]<-c('prop')
colnames(DiBS)[4]<-c('prop')
colnames(JLS)[4]<-c('prop')
colnames(LDS)[4]<-c('prop')
colnames(HBH)[4]<-c('prop')

TS$region<-'KT'
DiBS$region<-'KT'
JLS$region<-'KT'
LDS$region<-'KT'
HBH$region<-'KT'

GW$region<-'LD'
CK$region<-'LD'
SL$region<-'LD'
DaBS$region<-'LD'
GGB$region<-'LD'


CCG$region<-'LY'
YY$region<-'LY'
HR$region<-'LY'
TL$region<-'LY'
DC$region<-'LY'

JH$region<-'EC'
JMZ$region<-'EC'
HS$region<-'EC'
JQ$region<-'EC'
FNL$region<-'EC'

LD$region<-'NC'
BT$region<-'NC'
CJ$region<-'NC'
MYS$region<-'NC'
ZM$region<-'NC'

#Get longer table
composition_long<-rbind(TS, LDS, JLS, HBH, DiBS, 
     GW,CK, SL, DaBS, GGB,
     CCG, YY, HR, TL, DC,
     JH, JMZ, HS, JQ, FNL,
     LD, BT, CJ, MYS, ZM)
composition_long<-composition_long[c(2,4,5,6)]

#Get wider table
library(tidyr)
composition<-pivot_wider(composition_long, names_from = Mfgroup,
            values_from = prop)
composition[is.na(composition)]<-0
composition<-as.data.frame(composition)

#Further simplify the categories
composition$other<-composition$other+composition$others+composition$crioth #combine crinoids with other invertebrates
composition$others<-NULL
composition$crioth<-NULL
composition$cru<-composition$cur+composition$cru #cur is a wrong category, it needs to be combined with cru
composition$cur<-NULL 
composition$tws<-composition$wantws+composition$tws #combine all types of TWS together
composition$wantws<-NULL
#write.csv(composition, 'analysis_tables/Composition_AllSites.csv')

composition$cru_or_spoenc<-composition$cru+composition$spoenc #Since the difficulty to distinguish encrusting sponge and CCA, these two categories are merged as well.
composition$cru<-NULL
composition$spoenc<-NULL
write.csv(composition,'analysis_tables/Composition__CombineSpoencAndCru_AllSites.csv' )
#####Importing substrate data########
TS2<-read.csv('~/myExperiment/Kenting/TS_substrates_sum.csv')
LDS2<-read.csv('~/myExperiment/Kenting/LDS_substrates_sum.csv')
JLS2<-read.csv('~/myExperiment/Kenting/JLS_substrates_sum.csv')
HBH2<-read.csv('~/myExperiment/Kenting/HBH_substrates_sum.csv')
DiBS2<-read.csv('~/myExperiment/Kenting/DiBS_substrates_sum.csv')
GW2<-read.csv('~/myExperiment/Ludao/GW_substrates_sum.csv')
CK2<-read.csv('~/myExperiment/Ludao/CK_substrates_sum.csv')
SL2<-read.csv('~/myExperiment/Ludao/SL_substrates_sum.csv')
DaBS2<-read.csv('~/myExperiment/Ludao/DaBS_substrates_sum.csv')
GGB2<-read.csv('~/myExperiment/Ludao/GGB_substrates_sum.csv')

CCG2<-read.csv('~/myExperiment/Lanyu/CCG_substrates_sum.csv')
YY2<-read.csv('~/myExperiment/Lanyu/YY_substrates_sum.csv')
HR2<-read.csv('~/myExperiment/Lanyu/HR_substrates_sum.csv')
TL2<-read.csv('~/myExperiment/Lanyu/TL_substrates_sum.csv')
DC2<-read.csv('~/myExperiment/Lanyu/DC_substrates_sum.csv')

JH2<-read.csv('~/myExperiment/East Coast/JH_substrates_sum.csv')
JMZ2<-read.csv('~/myExperiment/East Coast/JMZ_substrates_sum.csv')
HS2<-read.csv('~/myExperiment/East Coast/HS_substrates_sum.csv')
JQ2<-read.csv('~/myExperiment/East Coast/JQ_substrates_sum.csv')
FNL2<-read.csv('~/myExperiment/East Coast/FNL_substrates_sum.csv')

LD2<-read.csv('~/myExperiment/North Coast/LD_substrates_sum.csv')
BT2<-read.csv('~/myExperiment/North Coast/BT_substrates_sum.csv')
CJ2<-read.csv('~/myExperiment/North Coast/CJ_substrates_sum.csv')
MYS2<-read.csv('~/myExperiment/North Coast/MYS_substrates_sum.csv')
ZM2<-read.csv('~/myExperiment/North Coast/82.5_substrates_sum.csv')


#Make longer table
substrates_long<-rbind(TS2, LDS2, JLS2, HBH2, DiBS2, 
      GW2,CK2, SL2, DaBS2, GGB2,
      CCG2, YY2, HR2, TL2, DC2,
      JH2, JMZ2, HS2, JQ2, FNL2,
      LD2, BT2, CJ2, MYS2, ZM2)
substrates_long<-substrates_long[c(2,4,5)]

#Make wider table
substrates<-pivot_wider(substrates_long, names_from = Mfgroup,
            values_from = prop)
substrates[is.na(substrates)]<-0
substrates<-as.data.frame(substrates)
substrates$region<-composition$region
substrates<-substrates[,-2]

#Save the preliminary composition table
write.csv(substrates,'analysis_tables/Substrates_AllSites.csv' )

3.1. Metric extraction (Fractal dimension)

This code is used to compute fractal dimension of different scale intervals based on surface area and planar area extracted from QGIS.

rm(list=ls())
setwd("~your working directory")
SA.data<-data.frame()
resolution<-c('The whole resolution', 1,2,4,8,16,32,64,128)
#Copy and paste surface area of different resolutions (accessed by function r.surf.area in GDAL)
SurfaceArea<-c( 'SA (whole Res)', 'SA (1 cm)', 'SA (2 cm)',
               'SA (4 cm)',   'SA (8 cm)',
                'SA (16 cm)',   'SA (32 cm)',
                'SA (64 cm)',   'SA (128 cm)')
#Copy and paste calculated planar area of different resolutions.
PlanAreaCalc<-c( 'PAc (Whole Res)', 
                 'PAc (1 cm)', 'PAc (2 cm)',
                 'PAc (4 cm)', 'PAc (8 cm)',
                 'PAc (16 cm)',  'PAc (32 cm)',
                 'PAc (64 cm)', ' PAc (128 cm)')
#Copy and paste null area (the true planar area is computed by subtracting calculated planar area from null area)
NullArea<-c( 'NA (Whole Res)', 
                 'NA (1 cm)', 'NA (2 cm)',
                 'NA (4 cm)', 'NA (8 cm)',
                 'NA (16 cm)',  'NA (32 cm)',
                 'NA (64 cm)', ' NA (128 cm)')
PlanArea<-PlanAreaCalc-NullArea #computing the true planar area
SA.data<-data.frame(resolution=resolution, 
           SurfaceArea=SurfaceArea, 
           PlanArea=PlanArea)
SA.data$SurfaceComplexity<-
  SA.data$SurfaceArea/SA.data$PlanArea
SA.data$LogSA<-log(SurfaceArea)
SA.data$LogRES<-log(resolution/100)
lmSA64<-lm(LogSA[2:8]~LogRES[2:8], data=SA.data) 
FD64<-2-lmSA64_2$coefficients[2] # Get D (1 - 64 cm)
library(ggplot2)
#Visualization of D (1 - 64 cm)
FD64Plot<-ggplot(SA.data[1:8, ], aes(x=LogRES, 
                                     y=LogSA))+
  geom_point()+
  geom_smooth(method='lm')+
  xlab('log(Resolution)')+
  ylab('log(Surface Area)')+
  ggtitle(paste('name of sampling site', 'D64= ', FD64_2))

####################################
#Get FD at different scale intervals
slopes<-c()
FD_intervals<-c()
for(i in 2:7){
  slopes[i-1]<-(SA.data$LogSA[i+1]-SA.data$LogSA[i])/(SA.data$LogRES[i+1]-SA.data$LogRES[i])
}
FD_intervals= 2-slopes
intervals<-c('D1_2','D2_4', 'D4_8', 'D8_16', 'D16_32', 'D32_64')
FDs<-data.frame(interv= intervals, FDs= FD_intervals)

3.2. Import metric table and conducting stepwise VIF selection on complexity metrics.

The metrics were collected by using QGIS. Fractal dimension (D) was computed by data extracted by QGIS and the code provided at section 1. All metrics of each site were manually organized in a table.

rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment")
library(corrplot)
library(car)
library(readxl)
library(olsrr)

#Read metric data
StrMetrics<-read.csv('Sites_MetricSummary.csv', header=T)
row.names(StrMetrics)<-StrMetrics[,1]
metrics<-StrMetrics[, 2:24]


#Import relief data
relief_data<-read_xlsx('VRandVRSD.xlsx', sheet=2)
relief_data<-as.data.frame(relief_data)
row.names(relief_data)<-relief_data[,1]
relief_data<-relief_data[-1]
metrics<-merge(metrics, relief_data[2], by='row.names')
row.names(metrics)<-metrics$Row.names
metrics<-metrics[-1] #Remove VR, keep VRSD (Sq)
write.csv(metrics, 'Site_MetricSummary_2.csv') #Save the table with all metrics
head(metrics)
##            S4         T4         V4      S32       T32        V32    PROC4
## 82.5 35.79695 0.10297550 0.03972595 34.08894 0.6829777 0.04316221 3.138004
## BT   29.06211 0.07379687 0.06354869 18.46841 0.3176737 0.01612316 4.993615
## CCG  31.98746 0.09315995 0.08002010 22.27737 0.3669287 0.01298073 4.518048
## CJ   37.43711 0.13200221 0.07057097 29.69820 0.6569097 0.08023318 3.669758
## CK   39.58948 0.13459857 0.05272271 28.06464 0.6684020 0.09568557 2.722089
## DaBS 22.30837 0.05852110 0.04697374 10.64565 0.2187278 0.01104392 3.920315
##           PLC4    PROC32    PLC32      S16       T16        V16    PROC16
## 82.5  9.242265 0.3542356 1.360989 34.22099 0.3514723 0.03748765 0.4714557
## BT   18.462807 0.3869983 1.543360 20.79732 0.2175307 0.03928283 1.0199474
## CCG  14.839175 0.2802628 1.027368 23.48149 0.2233868 0.02969492 0.8610387
## CJ   12.282814 0.4881684 1.838733 33.84235 0.3906941 0.07143632 0.9621377
## CK    8.519395 0.6309607 1.913845 34.12422 0.4080009 0.09667555 1.1375918
## DaBS 15.539315 0.2659814 1.967384 13.80435 0.1368370 0.01919955 0.7120210
##         PLC16      D64       SC     D1_2     D2_4     D4_8    D8_16   D16_32
## 82.5 5.380102 2.082288 1.637182 2.059848 2.059313 2.069394 2.048714 2.075094
## BT   5.696730 2.142668 1.617211 2.105806 2.136347 2.147361 2.146462 2.134220
## CCG  3.184027 2.141438 1.848779 2.196623 2.187020 2.149519 2.157569 2.034047
## CJ   3.296366 2.134910 1.979722 2.105106 2.102707 2.125147 2.165645 2.097948
## CK   3.775512 2.170443 2.080996 2.125388 2.149982 2.158457 2.136556 2.170571
## DaBS 3.822742 2.120286 1.374338 2.102362 2.106039 2.112706 2.097096 2.125669
##        D32_64     VRSD4
## 82.5 2.247949 0.6966732
## BT   2.187176 0.3594060
## CCG  2.140841 0.5719821
## CJ   2.238046 0.4688934
## CK   2.341126 0.3662118
## DaBS 2.214526 0.1714850
metrics_cor<-cor(metrics) #Get correlation matrix composed of every metric pair.

####################Corrplot adjustment####################################
metrics_pub<-metrics
colnames(metrics_pub)<-c('S - 4 cm', 'TRI - 4 cm', 'VRM - 4cm', 'S - 32 cm',
                         'TRI - 32 cm', 'VRM - 32 cm', 'PROC - 4 cm', 'PLC - 4 cm',
                         'PROC - 32 cm', 'PLC - 32 cm', 'S - 16 cm', 'TRI - 16 cm',
                         'VRM - 16 cm', 'PROC - 16 cm', 'PLC - 16 cm', 'D (1 - 64 cm)',
                         'SC', 'D (1 - 2 cm)', 'D (2 - 4 cm)', 'D (4 - 8 cm)', 'D (8 - 16 cm)',
                         'D (16 - 32 cm)', 'D (32 - 64 cm)', 'Sq')
metrics_cor<-cor(metrics_pub)
#corrplot(metrics_cor, method='number',
#         type='lower')
res<-cor.mtest(metrics_cor, 
               conf.level = .95)
#########################Displaying Corrplot#############################

corrplot(metrics_cor, method='number',
         type='lower', p.mat = res$p,
         sig.level = .05, tl.cex=0.5, cl.cex=0.5,
         number.cex=0.5)

#######################################################
###############Use usdm for stepwise VIF selection on metrics##################################
library(usdm)
v1<-vifstep(metrics, th=10)
metrics_selected<-exclude(metrics, v1)
metrics_selected #See metrics retained from selection
##              V4      PLC4    PROC32     PLC32    PLC16     D1_2     D2_4
## 82.5 0.03972595  9.242265 0.3542356 1.3609894 5.380102 2.059848 2.059313
## BT   0.06354869 18.462807 0.3869983 1.5433599 5.696730 2.105806 2.136347
## CCG  0.08002010 14.839175 0.2802628 1.0273676 3.184027 2.196623 2.187020
## CJ   0.07057097 12.282814 0.4881684 1.8387325 3.296366 2.105106 2.102707
## CK   0.05272271  8.519395 0.6309607 1.9138448 3.775512 2.125388 2.149982
## DaBS 0.04697374 15.539315 0.2659814 1.9673836 3.822742 2.102362 2.106039
## DC   0.07477847 17.502046 0.4166059 2.3384137 5.065466 2.175885 2.184300
## DiBS 0.04879224 13.948038 0.3448847 1.5995921 3.319096 2.146635 2.144978
## FNL  0.09231118 14.027762 0.5720152 1.8826472 3.643584 2.173418 2.185573
## GGB  0.06445081 14.743830 0.4614089 2.2478442 4.286704 2.142029 2.141366
## GW   0.04763611 15.339935 0.2874124 3.1911692 4.278659 2.125547 2.134165
## HBH  0.06507814 14.385386 0.5483947 1.7884567 3.607074 2.179179 2.169388
## HR   0.09157567 14.738506 0.5034706 2.2121285 4.682016 2.186643 2.206191
## HS   0.06418384 11.639053 0.7444362 2.4561036 3.899267 2.116186 2.117397
## JH   0.06854386 12.926524 0.5524131 1.9177405 3.275951 2.123935 2.120539
## JLS  0.08996626 13.711100 0.5674756 2.0121217 3.471119 2.213680 2.187870
## JMZ  0.10180321 13.580160 0.6749426 2.6838572 4.000385 2.147930 2.145423
## JQ   0.05014038 12.100439 0.5099651 2.0976212 2.811581 2.123612 2.102982
## LD   0.09773878 15.983127 0.2414503 0.8735963 4.529776 2.114213 2.142355
## LDS  0.07129418 18.708382 0.4316027 2.6060655 3.972630 2.117886 2.169715
## MYS  0.03591826 21.839115 0.1743400 1.7492916 4.682573 2.108662 2.116179
## SL   0.08072169 10.510869 0.5701214 1.8626170 3.517852 2.255038 2.152829
## TL   0.04717792 14.914221 0.2674686 1.8916910 4.247381 2.136491 2.124769
## TS   0.08013721 14.234224 0.4679286 2.5597380 4.639967 2.134288 2.130212
## YY   0.05256534 17.972904 0.2521757 1.9701828 5.083402 2.136654 2.139941
##         D8_16   D16_32   D32_64      VRSD4
## 82.5 2.048714 2.075094 2.247949 0.69667321
## BT   2.146462 2.134220 2.187176 0.35940599
## CCG  2.157569 2.034047 2.140841 0.57198211
## CJ   2.165645 2.097948 2.238046 0.46889336
## CK   2.136556 2.170571 2.341126 0.36621178
## DaBS 2.097096 2.125669 2.214526 0.17148500
## DC   2.144677 2.146597 2.195416 0.21552578
## DiBS 2.093201 2.055113 2.128465 0.27035985
## FNL  2.164326 2.171833 2.256766 0.32194592
## GGB  2.153132 2.114266 2.312947 0.35606082
## GW   2.092675 2.130449 2.179806 0.09572413
## HBH  2.094460 2.137286 2.292320 0.47868320
## HR   2.160204 2.106320 2.213140 0.29349338
## HS   2.185740 2.188789 2.338200 0.31157591
## JH   2.148578 2.156038 2.229615 0.49583007
## JLS  2.171962 2.132067 2.255832 0.51166963
## JMZ  2.197287 2.243912 2.308156 0.32193789
## JQ   2.138418 2.138883 2.238538 0.47133979
## LD   2.182984 2.173301 2.164168 0.49713277
## LDS  2.149087 2.092245 2.201565 0.17018665
## MYS  2.059969 2.003714 2.069643 0.10205495
## SL   2.106764 2.121427 2.303921 0.35146390
## TL   2.065154 2.127289 2.087826 0.26951771
## TS   2.128885 2.163362 2.132155 0.22094171
## YY   2.055983 2.084683 2.153331 0.17385214

4.1. Comparison of complexity metrics (Latitudinal correlation)

##Latitudinal correlation of metrics###
library(readxl)
#Import latitudes of sites
coordi<-read_xlsx('Photogrammetric_Plot_Coordinates.xlsx') 
## New names:
## • `` -> `...6`
coordi<-as.data.frame(coordi)
coordi<-coordi[,-c(4,5,6)]
row.names(coordi)<-coordi[,1]
#Merge latitudes and metrics
metrics_coordi<-merge(metrics, coordi,
                      by='row.names')
row.names(metrics_coordi)<-metrics_coordi[,1]
metrics_coordi<-metrics_coordi[-1]
metrics_coordi<-merge(metrics_coordi, StrMetrics[25],
                      by='row.names')
metrics_coordi$region<-factor(metrics_coordi$region,
                              levels=c('KT','LY','LD',
                                       'EC','NC'),
                              ordered = T)
#V4
cor_v4_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$V4, method='pearson',
         conf.level = 0.95)

#Plot the relationship with latitudes
library(ggplot2)
plot_v4_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=V4))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  annotate('text', x=22.5, y=0.037, size=5,label=paste('p = ',round(cor_v4_lat$p.value, 3)))+
  annotate('text', x=22.5, y=0.043, size=5,label=paste('r = ',round(cor_v4_lat$estimate, 3)))+
  xlab('Latitudes')+ ylab('VRM - 4 cm')+
  guides(color = guide_legend(title = "Region"))+
  scale_color_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_v4_lat
## `geom_smooth()` using formula = 'y ~ x'

ggsave('Figures/plot_v4_lat.jpeg', plot_v4_lat,
      width=15, height=15, unit='cm')
## `geom_smooth()` using formula = 'y ~ x'
#PROC32
cor_PROC32_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$PROC32, method='pearson',
                      conf.level = 0.95)

#Plot the relationship with latitudes
plot_PROC32_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=PROC32))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  annotate('text', x=24.5, y=0.135, size=5, label=paste('p = ',round(cor_PROC32_lat$p.value, 3)))+
  annotate('text', x=24.5, y=0.175, size=5, label=paste('r = ',round(cor_PROC32_lat$estimate, 3)))+
  xlab('Latitudes')+ ylab('PROC - 32 cm')+
  guides(color = guide_legend(title = "Region"))+
  scale_color_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_PROC32_lat
## `geom_smooth()` using formula = 'y ~ x'

ggsave('Figures/plot_PROC32_lat.jpeg', plot_PROC32_lat,
       width=15, height=15, unit='cm')
## `geom_smooth()` using formula = 'y ~ x'
#PLC4
cor_PLC4_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$PLC4, method='pearson',
                      conf.level = 0.95)

#Plot the relationship
plot_PLC4_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=PLC4))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  annotate('text', x=22.5, y=3, size=5,label=paste('p = ',round(cor_PLC4_lat$p.value, 3)))+
  annotate('text', x=22.5, y=5, size=5, label=paste('r = ',round(cor_PLC4_lat$estimate, 3)))+
  xlab('Latitudes')+ ylab('PLC - 4 cm')+
  guides(color = guide_legend(title = "Region"))+
  scale_color_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_PLC4_lat
## `geom_smooth()` using formula = 'y ~ x'

ggsave('Figures/plot_PLC4_lat.jpeg', plot_PLC4_lat,
       width=15, height=15, unit='cm')
## `geom_smooth()` using formula = 'y ~ x'
#PLC16
cor_PLC16_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$PLC16, method='pearson',
                     conf.level = 0.95)

#Plot relationship
plot_PLC16_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=PLC16))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  annotate('text', x=22.5, y=2.9,size=5, label=paste('p = ',round(cor_PLC16_lat$p.value, 3)))+
  annotate('text', x=22.5, y=3.2, size=5,label=paste('r = ',round(cor_PLC16_lat$estimate, 3)))+
  xlab('Latitudes')+ ylab('PLC - 16 cm')+
  guides(color = guide_legend(title = "Region"))+
  scale_color_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_PLC16_lat
## `geom_smooth()` using formula = 'y ~ x'

ggsave('Figures/plot_PLC16_lat.jpeg', plot_PLC16_lat,
       width=15, height=15, unit='cm')
## `geom_smooth()` using formula = 'y ~ x'
#PLC32
cor_PLC32_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$PLC32, method='pearson',
                                              conf.level = 0.95)

#no significant correlation btw sc and latitude
plot_PLC32_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=PLC32))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  annotate('text', x=22.6, y=1, size=5,label=paste('p = ',round(cor_PLC32_lat$p.value, 3)))+
  annotate('text', x=22.6, y=1.2, size=5,label=paste('r = ',round(cor_PLC32_lat$estimate, 3)))+
  xlab('Latitudes')+ ylab('PLC - 32 cm')+
  guides(color = guide_legend(title = "Region"))+
  scale_color_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_PLC32_lat
## `geom_smooth()` using formula = 'y ~ x'

ggsave('figures/plot_PLC32_lat.jpeg', plot_PLC32_lat,
       width=15, height=15, unit='cm')
## `geom_smooth()` using formula = 'y ~ x'
#D12

cor_D12_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$D1_2, method='pearson',
                     conf.level = 0.95)

#Plot relationship
plot_D12_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=D1_2))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  annotate('text', x=22.5, y=2.07,size=6.5, label=paste('p = ',round(cor_D12_lat$p.value, 3)))+
  annotate('text', x=22.5, y=2.085, size=6.5,label=paste('r = ',round(cor_D12_lat$estimate, 3)))+
  xlab('Latitudes')+ ylab('D (1 - 2 cm)')+
  guides(color = guide_legend(title = "Region"))+
  scale_color_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_D12_lat
## `geom_smooth()` using formula = 'y ~ x'

ggsave('Figures/plot_D12_lat.jpeg', plot_D12_lat,
       width=15, height=15, unit='cm')
## `geom_smooth()` using formula = 'y ~ x'
#D24
cor_D24_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$D2_4, method='pearson',
                        conf.level = 0.95)

#Plot relationship
plot_D24_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=D2_4))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  annotate('text', x=22.5, y=2.07, size=6.5, label=paste('p = ',round(cor_D24_lat$p.value, 3)))+
  annotate('text', x=22.5, y=2.085, size=6.5, label=paste('r = ',round(cor_D24_lat$estimate, 3)))+
  xlab('Latitudes')+ ylab('D (2 - 4 cm)')+
  guides(color = guide_legend(title = "Region"))+
  scale_color_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_D24_lat
## `geom_smooth()` using formula = 'y ~ x'

ggsave('Figures/plot_D24_lat.jpeg', plot_D24_lat,
       width=15, height=15, unit='cm')
## `geom_smooth()` using formula = 'y ~ x'
#D816

cor_D816_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$D8_16, method='pearson',
                      conf.level = 0.95)

#Plot relationship
plot_D816_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=D8_16))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  annotate('text', x=23, y=2.057, size=5,label=paste('p = ',round(cor_D816_lat$p.value, 3)))+
  annotate('text', x=23, y=2.07, size=5,label=paste('r = ',round(cor_D816_lat$estimate, 3)))+
  xlab('Latitudes')+ ylab('D (8 - 16 cm)')+
  guides(color = guide_legend(title = "Region"))+
  scale_color_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_D816_lat
## `geom_smooth()` using formula = 'y ~ x'

ggsave('Figures/plot_D816_lat.jpeg', plot_D816_lat,
       width=15, height=15, unit='cm')
## `geom_smooth()` using formula = 'y ~ x'
#D1632
cor_D1632_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$D16_32, method='pearson',
                      conf.level = 0.95)


#Plot relationship
plot_D1632_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=D16_32))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  annotate('text', x=23, y=2.02,size=5, label=paste('p = ',round(cor_D1632_lat$p.value, 3)))+
  annotate('text', x=23, y=2.04, size=5,label=paste('r = ',round(cor_D1632_lat$estimate, 3)))+
  xlab('Latitudes')+ ylab('D (16 - 32 cm)')+
  guides(color = guide_legend(title = "Region"))+
  scale_color_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_D1632_lat
## `geom_smooth()` using formula = 'y ~ x'

ggsave('Figures/plot_D1632_lat.jpeg', plot_D1632_lat,
       width=15, height=15, unit='cm')
## `geom_smooth()` using formula = 'y ~ x'
#D3264
cor_D3264_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$D32_64, method='pearson',
                      conf.level = 0.95)

#Plot relationship
plot_D3264_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=D32_64))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  annotate('text', x=23, y=2.07, size=5,label=paste('p = ',round(cor_D3264_lat$p.value, 3)))+
  annotate('text', x=23, y=2.09, size=5,label=paste('r = ',round(cor_D3264_lat$estimate, 3)))+
  xlab('Latitudes')+ ylab('D (32 - 64 cm)')+
  guides(color = guide_legend(title = "Region"))+
  scale_color_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_D3264_lat
## `geom_smooth()` using formula = 'y ~ x'

ggsave('Figures/plot_D3264_lat.jpeg', plot_D3264_lat,
       width=15, height=15, unit='cm')
## `geom_smooth()` using formula = 'y ~ x'
#Sq

cor_VRSD_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$VRSD4, method='pearson',
                       conf.level = 0.95)

#Plot relationship
plot_VRSD_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=VRSD4))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  annotate('text', x=23.5, y=0.15, size=5,label=paste('p = ',round(cor_VRSD_lat$p.value, 3)))+
  annotate('text', x=23.5, y=0.18, size=5,label=paste('r = ',round(cor_VRSD_lat$estimate, 3)))+
  labs(y='Sq')+xlab('Latitudes')+
  guides(color = guide_legend(title = "Region"))+
  scale_color_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_VRSD_lat
## `geom_smooth()` using formula = 'y ~ x'

ggsave('Figures/plot_VRSD_lat.jpeg', plot_VRSD_lat,
       width=15, height=15, unit='cm')
## `geom_smooth()` using formula = 'y ~ x'

4.2. Comparison of complexity metrics (Interregional differences)

##########################Interregional comparisons: Kruskal Wallis test##############################


V4_kruskal<-kruskal.test(V4~region, data=metrics_coordi)
V32_kruskal<-kruskal.test(V32~region, data=metrics_coordi)
PROC32_kruskal<-kruskal.test(PROC32~region, data=metrics_coordi)
PLC4_kruskal<-kruskal.test(PLC4~region, data=metrics_coordi)
PLC16_kruskal<-kruskal.test(PLC16~region, data=metrics_coordi)
PLC32_kruskal<-kruskal.test(PLC32~region, data=metrics_coordi)
D12_kruskal<-kruskal.test(D1_2~region, data=metrics_coordi)
D24_kruskal<-kruskal.test(D2_4~region, data=metrics_coordi)
D48_kruskal<-kruskal.test(D4_8~region, data=metrics_coordi)
D816_kruskal<-kruskal.test(D8_16~region, data=metrics_coordi)
D1632_kruskal<-kruskal.test(D16_32~region, data=metrics_coordi)
D3264_kruskal<-kruskal.test(D32_64~region, data=metrics_coordi)
VRSD_kruskal<-kruskal.test(VRSD4~region, data=metrics_coordi)

######################Pairwise Dunn test########################################
library(PMCMRplus)
PROC32_DT<-kwAllPairsDunnTest(PROC32~region, data=metrics_coordi, method='bh')
#Significance at EC-LY, EC-NC
D3264_DT<-kwAllPairsDunnTest(D32_64~region, data=metrics_coordi, method="bh")
#Significance may exist at EC-LY
D12_DT<-kwAllPairsDunnTest(D1_2~region, data=metrics_coordi, method="bh")
#Significance at LY-NC, KT-NC

library(rcompanion)
library(rempsyc)
#Save table
PROC32_DT<-PMCMRTable(PROC32_DT)
D3264_DT<-PMCMRTable(D3264_DT)
D12_DT<-PMCMRTable(D12_DT)
DT_sum<-cbind(D12_DT, D3264_DT, PROC32_DT)
DT_sum[c(3,5)]<-NULL
colnames(DT_sum)[2:4]<-c('p.value of D (1 - 2 cm) cm', 'p.value of D (32 - 64 cm)', 'p.value of PROC - 32 cm')
DT_NT<-nice_table(DT_sum)
DT_NT

Comparison

p.value of D (1 - 2 cm) cm

p.value of D (32 - 64 cm)

p.value of PROC - 32 cm

LY - KT = 0

1

1

0.928

LD - KT = 0

1

0.72

1

EC - KT = 0

1

0.563

0.928

NC - KT = 0

0.0537

1

0.928

LD - LY = 0

1

0.163

0.928

EC - LY = 0

1

0.0994

0.0411

NC - LY = 0

0.0127

1

1

EC - LD = 0

1

1

0.75

NC - LD = 0

0.348

0.547

0.928

NC - EC = 0

0.411

0.385

0.0263

flextable::save_as_docx(DT_NT, path = "Figures/PairwiseComparison_indicators_Dunn.docx")


####Plot relationships between regions###
library(ggpubr)
#V4
box_V4_region<-ggplot(data=metrics_coordi, 
                      aes(x=region, y=V4 ))+
  geom_boxplot(aes(fill=region))+
 
  labs(y='VRM - 4 cm')+xlab('Region')+
  guides(fill = guide_legend(title = "Region"))+
  scale_fill_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))+
  stat_compare_means(label.y = 0.1, cex=5)
box_V4_region

ggsave('Figures/box_v4_region.jpeg', box_V4_region,
       width=15, height=15, unit='cm')

#PROC32

box_PROC32_region<-ggplot(data=metrics_coordi, 
                       aes(x=region, y=PROC32 ))+
  geom_boxplot(aes(fill=region))+
  ggtitle(label = paste('KW chi-squared= ', round(PROC32_kruskal$statistic, 3)))+
  labs(y='PROC - 32 cm')+xlab('Region')+
  guides(fill = guide_legend(title = "Region"))+
  scale_fill_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))


box_PROC32_region<-box_PROC32_region+
  stat_compare_means(label.y = 0.17, cex=5)+ # Add global p-value
  annotate('text', x=1:5, y=0.75, label=c('ns', '*', 'ns', '', '*'), 
           cex=5) 
box_PROC32_region

ggsave('Figures/box_PROC32_region.jpeg', box_PROC32_region,
       width=15, height=15, unit='cm')


#PLC4

box_PLC4_region<-ggplot(data=metrics_coordi, 
                         aes(x=region, y=PLC4 ))+
  geom_boxplot(aes(fill=region))+
  ggtitle(label =paste('KW chi-squared= ', round(PLC4_kruskal$statistic, 3)))+
  labs(y='PLC - 4 cm')+xlab('Region')+
  guides(fill = guide_legend(title = "Region"))+
  scale_fill_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))+
  stat_compare_means(label.y = 21, cex=5)
box_PLC4_region

ggsave('Figures/box_PLC4_region.jpeg', box_PLC4_region,
       width=15, height=15, unit='cm')


#PLC16
box_PLC16_region<-ggplot(data=metrics_coordi, 
                      aes(x=region, y=PLC16 ))+
  geom_boxplot(aes(fill=region))+
  ggtitle(label =paste('KW chi-squared= ', round(PLC16_kruskal$statistic, 3)))+
  labs(y='PLC - 16 cm')+xlab('Region')+
  guides(fill = guide_legend(title = "Region"))+
  scale_fill_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))+
  stat_compare_means(label.y = 5.5, cex=5)

box_PLC16_region

ggsave('Figures/box_PLC16_region.jpeg', box_PLC16_region,
       width=15, height=15, unit='cm')
#PLC32
box_PLC32_region<-ggplot(data=metrics_coordi, 
                      aes(x=region, y=PLC32 ))+
  geom_boxplot(aes(fill=region))+
  ggtitle(label =paste('KW chi-squared= ', round(PLC32_kruskal$statistic, 3)))+
  labs(y='PLC - 32 cm')+xlab('Region')+
  guides(fill = guide_legend(title = "Region"))+
  scale_fill_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))+
  stat_compare_means(label.y =3.15, cex=4)

box_PLC32_region

ggsave('Figures/box_PLC32_region.jpeg', box_PLC32_region,
       width=15, height=15, unit='cm')
#D12
box_D12_region<-ggplot(data=metrics_coordi, 
                      aes(x=region, y=D1_2 ))+
  geom_boxplot(aes(fill=region))+
  ggtitle(label =paste('KW chi-squared= ', round(D12_kruskal$statistic, 3)))+
  labs(y='D(1 - 2 cm)')+xlab('Region')+
  guides(fill = guide_legend(title = "Region"))+
  scale_fill_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))


box_D12_region<-box_D12_region+
  stat_compare_means(label.y = 2.32, cex=5)+ # Add global p-value
  annotate('text', x=1:5, y=2.28, label=c('ns', '*', 'ns','ns',''), 
           cex=5)

box_D12_region

ggsave('Figures/box_D12_region.jpeg', box_D12_region,
       width=15, height=15, unit='cm')


#D24

box_D24_region<-ggplot(data=metrics_coordi, 
                       aes(x=region, y=D2_4 ))+
  geom_boxplot(aes(fill=region))+
  ggtitle(label =paste('KW chi-squared= ', round(D24_kruskal$statistic, 3)))+
  labs(y='D(2 - 4 cm) cm')+xlab('Region')+
  guides(fill = guide_legend(title = "Region"))+
  scale_fill_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))+
  stat_compare_means(label.y = 2.06, cex=5)

box_D24_region

ggsave('Figures/box_D24_region.jpeg', box_D24_region,
       width=15, height=15, unit='cm')


#D816
box_D816_region<-ggplot(data=metrics_coordi, 
                       aes(x=region, y=D8_16 ))+
  geom_boxplot(aes(fill=region))+
  ggtitle(label =paste('KW chi-squared= ', round(PLC4_kruskal$statistic, 3)))+
  labs(y='D(8 - 16 cm)')+xlab('Region')+
  guides(fill = guide_legend(title = "Region"))+
  scale_fill_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))+
  stat_compare_means(label.y = 2.19, cex=5)

box_D816_region

ggsave('Figures/box_D816_region.jpeg', box_D816_region,
       width=15, height=15, unit='cm')
#D1632
box_D1632_region<-ggplot(data=metrics_coordi, 
                        aes(x=region, y=D16_32 ))+
  geom_boxplot(aes(fill=region))+
  ggtitle(label =paste('KW chi-squared= ', round(PLC4_kruskal$statistic, 3)))+
  labs(y='D(16 - 32 cm)')+xlab('Region')+
  guides(fill = guide_legend(title = "Region"))+
  scale_fill_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))+
  stat_compare_means(label.y = 2.23, cex=5)

box_D1632_region

ggsave('Figures/box_D1632_region.jpeg', box_D1632_region,
       width=15, height=15, unit='cm')
#D3264
box_D3264_region<-ggplot(data=metrics_coordi, 
                         aes(x=region, y=D32_64 ))+
  geom_boxplot(aes(fill=region))+
  ggtitle(label =paste('KW chi-squared= ', round(PLC4_kruskal$statistic, 3)))+
  labs(y='D(32 - 64 cm)')+xlab('Region')+
  guides(fill = guide_legend(title = "Region"))+
  scale_fill_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))

box_D3264_region<-box_D3264_region+
  stat_compare_means(label.y = 2.4, cex=5)+ # Add global p-value
  # Pairwise comparison against reference
  annotate('text', x=1:5, y=2.36, label=c('ns', 'ns', 'ns', '', 'ns'),
           cex=5)
box_D3264_region

ggsave('Figures/box_D3264_region.jpeg', box_D3264_region,
       width=15, height=15, unit='cm')
#Sq
box_VRSD_region<-ggplot(data=metrics_coordi, 
                         aes(x=region, y=VRSD4 ))+
  geom_boxplot(aes(fill=region))+
  ggtitle(label =paste('KW chi-squared= ', round(PLC4_kruskal$statistic, 3)))+
  labs(y='Sq')+xlab('Region')+
  guides(fill = guide_legend(title = "Region"))+
  scale_fill_manual(values=c( "#fc0f00",  "#fb7200","#fec700", "#15e0fa","#0c59fe" ))

box_VRSD_region<-box_VRSD_region+
  stat_compare_means(label.y = 0.65, cex=5)
  
box_VRSD_region

ggsave('Figures/box_VRSD_region.jpeg', box_VRSD_region,
       width=15, height=15, unit='cm')

4.3. PCA based on complexity metrics

PCA was used to visualize the relationships between sites, PERMANOVA was used to examined the interregional differences of complexity.

###################PCA based on metrics#########################
library(vegan)
## 載入需要的套件:permute
## 載入需要的套件:lattice
## This is vegan 2.6-4
metr_sel_std<-decostand(metrics_selected, 'standardize')
metrics_pca<-rda(metr_sel_std, scale=T)
PCAsummary<-summary(metrics_pca)

gps<-c('NC', 'NC', 'LY', 'NC', 'LD', 'LD', 'LY', 'KT', 'EC', 'LD', 'LD', 'KT', 'LY',
       'EC', 'EC', 'KT', 'EC', 'EC', 'NC', 'KT', 'NC', 'LD', 'LY', 'KT', 'LY')
my_cols<-c("#15e0fa", "#fc0f00","#fec700",  "#fb7200" ,"#0c59fe" )
gps2<-as.factor(gps)
levels(gps2)<-my_cols
#Visualization
ordiplot(metrics_pca, type='n',display='site',
         xlab=paste('PC1,', 
                    round(PCAsummary[["cont"]][["importance"]][3,1]*100, 2), '%'),
         ylab=paste('PC2,', 
                    round(PCAsummary[["cont"]][["importance"]][2,2]*100, 2), '%'),
         scaling=2)
pca.scores<-scores(metrics_pca, 
                   choices=1:2, 
                   scaling=2,
                   display='site')
points(x=pca.scores[,1], y=pca.scores[,2], pch=20, cex=1.5,
       col=as.vector(gps2))

ordispider(metrics_pca, 
           groups=gps,
           display='site',
           label=F,
           lwd=2,
           col=c("#15e0fa", "#fc0f00","#fec700",  "#fb7200" ,"#0c59fe" ),
           scaling=2)

ordihull(metrics_pca,
         groups=gps,
         display='site',
         draw='lines',
        lwd = 2,
         col=c("#15e0fa", "#fc0f00","#fec700",  "#fb7200" ,"#0c59fe" ),
         scaling=2)
vct.scores<-scores(metrics_pca, 
                   choices=1:2, 
                   scaling=2,
                   display='sp')
arrows(x0=0, y0=0,
       x1=vct.scores[,1],
       y1=vct.scores[,2],
       length=0.1,
       lty=1,
       lwd=3,
       col=rgb(0.8, 0.1,0,0.7))
text(x=vct.scores[1,1]+0.3, y= vct.scores[1,2]+0.2,
'VRM - 4 cm', col= 'red', cex=1)
text(x=vct.scores[2,1]-0.3, y= vct.scores[2,2]+0.2,
     'PLC - 4 cm', col= 'red', cex=1)
text(x=vct.scores[3,1]+0.5, y= vct.scores[3,2],
     'PROC - 32 cm', col= 'red', cex=1)
text(x=vct.scores[4,1]+0.28, y= vct.scores[4,2]+0.26,
     'PLC - 32 cm', col= 'red', cex=0.8)
text(x=vct.scores[5,1]-0.8, y= vct.scores[5,2]+0.4,
     'PLC - 16 cm', col= 'red', cex=1)
text(x=vct.scores[6,1]+0.3, y= vct.scores[6,2]+0.3,
     'D (1 - 2 cm)', col= 'red', cex=1)
text(x=vct.scores[7,1]+0.3, y= vct.scores[7,2]+0.2,
     'D (2 - 4 cm)', col= 'red', cex=1)
text(x=vct.scores[8,1]+0.62, y= vct.scores[8,2]+0.1,
     'D (8 - 16 cm)', col= 'red', cex=1)
text(x=vct.scores[9,1]+0.75, y= vct.scores[9,2],
     'D (16 - 32 cm)', col= 'red', cex=1)
text(x=vct.scores[10,1]+0.55, y= vct.scores[10,2]-0.1,
     'D (32 - 64 cm)', col= 'red', cex=1)
text(x=vct.scores[11,1]+0.2, y= vct.scores[11,2]-0.4,
     'Sq', col= 'red', cex=1)
legend('topleft', fill=my_cols,
       legend= c('East Coast','Kenting', 'Ludao', 'Lanyu', 'North Coast'), 
       title='Region', cex=1)

#PERMANOVA on complexity
metr_sel_euc<-vegdist(metr_sel_std, 'euclidean')

metr_sel_region<-cbind(metr_sel_std, gps)

metr_sel_PERMANOVA<-adonis2(metr_sel_std~gps, data=metr_sel_region,
                            permutations = 999, method='euclidean')

metr_sel_PERMANOVA
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = metr_sel_std ~ gps, data = metr_sel_region, permutations = 999, method = "euclidean")
##          Df SumOfSqs      R2      F Pr(>F)   
## gps       4   79.403 0.30077 2.1507  0.007 **
## Residual 20  184.597 0.69923                 
## Total    24  264.000 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
PERMANOVA_nice<-nice_table(metr_sel_PERMANOVA)
flextable::save_as_docx(PERMANOVA_nice, path = "Figures/Complexity_PERMANOVA.docx")
library(RVAideMemoire)
## *** Package RVAideMemoire v 0.9-81-2 ***
## 
## 載入套件:'RVAideMemoire'
## 下列物件被遮斷自 'package:raster':
## 
##     cv
library(devtools)
## 載入需要的套件:usethis
## 
## 載入套件:'devtools'
## 下列物件被遮斷自 'package:permute':
## 
##     check
library(pairwiseAdonis)
## 載入需要的套件:cluster
#Pairwise permutation test
permtest<-pairwise.adonis(as.data.frame(metr_sel_std), metr_sel_region$gps,
                          sim.method='euclidean')

permtest_nice<-nice_table(permtest,note='* p < .05, ** p < .01, *** p < .001', highlight = T)
permtest_nice

pairs

Df

SumsOfSqs

F.Model

R2

p

p.adjusted

sig

NC vs LY

1.00

18.07

1.49

.16

.179

1.00

NC vs LD

1.00

24.46

2.04

.20

.050

0.50

NC vs KT

1.00

22.74

2.01

.20

.076

0.76

NC vs EC

1.00

40.54

3.90

.33

.008**

0.08

LY vs LD

1.00

16.27

1.82

.19

.120

1.00

LY vs KT

1.00

6.09

0.74

.08

.590

1.00

LY vs EC

1.00

35.73

4.84

.38

.008**

0.08

LD vs KT

1.00

6.85

0.84

.10

.485

1.00

LD vs EC

1.00

13.04

1.81

.18

.110

1.00

KT vs EC

1.00

14.72

2.25

.22

.067

0.67

Note. * p < .05, ** p < .01, *** p < .001

flextable::save_as_docx(permtest_nice, path = "Figures/Complexity_PairwisePermTest.docx")

5. Further simplification on composition data and conducting PCA

Several categories were further merged (i.e., combine Heliopora, Millepora, Tubipora with hard corals). And combining biotic and abiotic (substrates) datasets again. Using PCA to visualize relative relationships between sites, using PERMANOVA to examine the interregional differences.

rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment")
library(corrplot)
library(car)
library(readxl)
library(olsrr)

#Read metric data
StrMetrics<-read.csv('Sites_MetricSummary.csv', header=T)
row.names(StrMetrics)<-StrMetrics[,1]
metrics<-StrMetrics[, 2:24]


#Import relief data
relief_data<-read_xlsx('VRandVRSD.xlsx', sheet=2)
relief_data<-as.data.frame(relief_data)
row.names(relief_data)<-relief_data[,1]
relief_data<-relief_data[-1]
metrics<-merge(metrics, relief_data[2], by='row.names')
row.names(metrics)<-metrics$Row.names
metrics<-metrics[-1]
write.csv(metrics, 'Site_MetricSummary_2.csv')

#########################################################################
########Import composition (Both biotic components and substrates) and further organization###########################
library(dplyr)
## 
## 載入套件:'dplyr'
## 下列物件被遮斷自 'package:raster':
## 
##     intersect, select, union
## 下列物件被遮斷自 'package:car':
## 
##     recode
## 下列物件被遮斷自 'package:stats':
## 
##     filter, lag
## 下列物件被遮斷自 'package:base':
## 
##     intersect, setdiff, setequal, union
composition<-read.csv('analysis_tables/Composition__CombineSpoencAndCru_AllSites.csv',row.names = 1 ) #Biotic components
substrates<-read.csv( 'analysis_tables/Substrates_AllSites.csv') #substrates
substrates<-substrates[,-1]
composition_substrates<-cbind(composition, substrates[,-c(1,6)])
composition_substrates$bss<-1-rowSums(composition_substrates[,3:45]) #get BSS
composition_substrates$hcbus<-composition_substrates$hcbus+composition_substrates$helbus+
  composition_substrates$milbus+composition_substrates$hcdig #Get bushy HC

composition_substrates$helbus<-NULL
composition_substrates$milbus<-NULL
composition_substrates$hcdig<-NULL

#get encrusting HC
composition_substrates$hcenc<-composition_substrates$hcenc+
  composition_substrates$milenc+composition_substrates$helenc+
  composition_substrates$tubenc
composition_substrates$milenc<-NULL
composition_substrates$helenc<-NULL
composition_substrates$tubenc<-NULL

#get massive HC
composition_substrates$hcmas<-composition_substrates$hcmas+composition_substrates$tubmas
composition_substrates$tubmas<-NULL

#get Arborescent HC
composition_substrates$hcarb<-composition_substrates$hcarb+composition_substrates$milarb
composition_substrates$milarb<-NULL

#get macroalgae
composition_substrates$ma<-composition_substrates$com+
  composition_substrates$art+composition_substrates$cof
composition_substrates$com<-NULL
composition_substrates$art<-NULL
composition_substrates$cof<-NULL

#get US
composition_substrates$us<-composition_substrates$sandus+
  composition_substrates$rubbus+composition_substrates$siltus+
  composition_substrates$gravus
composition_substrates$sandus<-NULL
composition_substrates$rubbus<-NULL
composition_substrates$siltus<-NULL
composition_substrates$gravus<-NULL

head(composition_substrates)
##   site region      actinia      ceralg        hcarb      hcbus       hccol
## 1   TS     KT 1.382217e-03 0.001840512 7.365073e-03 0.07568280 0.001178526
## 2  LDS     KT 0.000000e+00 0.000000000 1.601930e-02 0.11848354 0.014356737
## 3  JLS     KT 4.594236e-05 0.000148033 5.171541e-02 0.13114650 0.000000000
## 4  HBH     KT 0.000000e+00 0.000000000 0.000000e+00 0.05456729 0.007572080
## 5 DiBS     KT 0.000000e+00 0.000000000 5.505862e-05 0.07930534 0.002333195
## 6   GW     LD 0.000000e+00 0.000000000 0.000000e+00 0.06344037 0.000000000
##        hcenc        hcfol      hcmas       hctab        ocdig        oclob
## 1 0.16715161 0.0050774656 0.06404287 0.040455649 0.0001649510 3.451637e-03
## 2 0.31648276 0.0313315459 0.04485327 0.000000000 0.0124688964 2.069143e-03
## 3 0.03559424 0.0423277998 0.00963379 0.002503439 0.0000000000 1.443082e-03
## 4 0.04454544 0.4555850206 0.01511464 0.000000000 0.0000000000 3.233330e-03
## 5 0.10917749 0.0002641661 0.01333520 0.000000000 0.0003398285 8.729073e-05
## 6 0.18696224 0.2107142169 0.02108280 0.000000000 0.0002013752 0.000000e+00
##          ocmas          tws          fil       ocbus    ascidian        ocfan
## 1 0.0021837048 3.937554e-04 0.0000000000 0.000000000 0.000000000 0.000000e+00
## 2 0.0009843722 0.000000e+00 0.0018783695 0.005894327 0.000000000 0.000000e+00
## 3 0.0000000000 0.000000e+00 0.0035133109 0.002369318 0.001000274 4.835259e-05
## 4 0.0000000000 0.000000e+00 0.0006558537 0.005005998 0.000000000 0.000000e+00
## 5 0.0002757979 9.899104e-05 0.0000000000 0.007079813 0.000000000 0.000000e+00
## 6 0.0000000000 3.013859e-04 0.0000000000 0.000000000 0.000000000 0.000000e+00
##         zoaenc       spoglo       spomas       zoamas     occlu ocenc
## 1 0.0000000000 0.000000e+00 0.0000000000 0.0000000000 0.0000000     0
## 2 0.0000000000 0.000000e+00 0.0000000000 0.0000000000 0.0000000     0
## 3 0.0002374461 0.000000e+00 0.0000000000 0.0000000000 0.0000000     0
## 4 0.0000000000 0.000000e+00 0.0000000000 0.0000000000 0.0000000     0
## 5 0.0008747468 2.266706e-05 0.0001298605 0.0002267812 0.0000000     0
## 6 0.0000000000 0.000000e+00 0.0000000000 0.0009696243 0.1732429     0
##          other sporep spopap coeenc cru_or_spoenc       bss           ma
## 1 3.117624e-03      0      0      0  0.0004121527 0.4972918 0.0681254510
## 2 8.293761e-05      0      0      0  0.0058026699 0.4038181 0.0254740217
## 3 7.965084e-05      0      0      0  0.0057743596 0.6127960 0.0830581394
## 4 0.000000e+00      0      0      0  0.0009811873 0.1511659 0.0010077042
## 5 1.610879e-03      0      0      0  0.0101201119 0.6875589 0.0693698048
## 6 0.000000e+00      0      0      0  0.0091574471 0.3329709 0.0009566711
##           us
## 1 0.06068224
## 2 0.00000000
## 3 0.01656496
## 4 0.26056556
## 5 0.01773404
## 6 0.00000000
##Combine metrics with composition
#Import metrics
metrics<-read.csv('Site_MetricSummary_2.csv', row.names = 1, header=T)
row.names(composition_substrates)<-composition_substrates$site
comp_sub_metrics<-merge(composition_substrates, metrics, by='row.names')
row.names(comp_sub_metrics)<-comp_sub_metrics$Row.names
comp_sub_metrics<-comp_sub_metrics[,-1]
#Change colnames
indicators<-c('S - 4 cm', 'TRI - 4 cm', 'VRM - 4 cm', 
              'S - 32 cm', 'TRI - 32 cm', 'VRM - 32 cm', 
              'PROC - 4 cm', 'PLC - 4 cm', 'PROC - 32 cm',
              'PLC - 32 cm', 'S - 16 cm', 'TRI - 16 cm', 'VRM - 16 cm',
              'PROC - 16 cm', 'PLC - 16 cm', 'D(1 - 64 cm)', 'SC',
              'D (1 - 2 cm)', 'D (2 - 4 cm)', 'D (4 - 8 cm)', 
              'D(8 - 16 cm)', 'D(16 - 32 cm)', 'D(32 - 64 cm)',
              'Sq')
colnames(comp_sub_metrics)[34:57]<-indicators


#####################PCA for composition#########################################
##############set color###########################

library(vegan)
library(fishualize)
library(devtools)
#devtools::install_github("nschiett/fishualize", force = TRUE)

my_cols<-c("#15e0fa", "#fc0f00","#fec700",  "#fb7200" ,"#0c59fe" )

##########################tb-PCA (hellinger transformed)###########################################
#Hellinger transformation on composition
comp_sub_hel<-decostand(composition_substrates[,3:33], method='hellinger')
comp_sub_hel<-decostand(comp_sub_hel, method='standardize')
colnames(comp_sub_hel)<-c('Actinia', 'Ceratodictyon', 'Arborescent HC',
                          'Bushy HC', 'Columnar HC', 'Encrusting HC',
                          'Foliose HC', 'Massive HC', 'Tabular HC', 'Digitate OC',
                          'Lobate OC', 'Massive OC', 'TWS','Filamentous algae', 'Bushy OC',
                          'Ascidian','Fan OC', 'Encrusting Zoanthids',
                          'Globular Sponges', 'Massive Sponges','Massive Zoanthids',
                          'Clustered OC','Encrusting OC', 'Others', 'Repent Sponges',
                          'Papillate Sponges', 'Corallimorpharia', 'CCA',
                          'BSS', 'Macroalgae', 'US')
groups<-composition_substrates$region
comp_sub_RDA<-rda(comp_sub_hel)
RDA_summary<-summary(comp_sub_RDA)
#Visualization of PCA

ordiplot(comp_sub_RDA, type='n',
         display=c('site', 'sp'),
         scaling=3,
         xlab=paste('PC1, ', 
                    round(RDA_summary$cont$importance[2,1]*100, 2), '%'),
         ylab=paste('PC2, ',round(RDA_summary$cont$importance[2,2]*100, 2), '%'),
        )
title(sub='Benthic composition based on sites (hellinger transformed)')
ordihull(comp_sub_RDA, 
         groups=groups,
         draw='lines', 
         scaling=3,
         col=my_cols,
         lwd=2,
         alpha=0.8)
vct.scores<-scores(comp_sub_RDA, 
                   choices=1:2, 
                   display='sp',
                   scaling=3)

arrows(x0=0, y0=0,
       x1=vct.scores[,1],
       y1=vct.scores[,2],
       length=0,
       lty=1,
       lwd=0.7,
       col='red')
text(vct.scores[3,1],vct.scores[3,2]-0.05,
     row.names(vct.scores)[3],
     col='red', cex=0.8) #Arborescent HC
text(vct.scores[4,1]+0.13,vct.scores[4,2],
     row.names(vct.scores)[4],
     col='red', cex=0.8) #Bushy HC
text(vct.scores[5,1]+0.11,vct.scores[5,2]-0.02,
     row.names(vct.scores)[5],
     col='red', cex=0.6) #Columnar HC
text(vct.scores[6,1]+0.03,vct.scores[6,2]-0.06,
     row.names(vct.scores)[6],
     col='red', cex=0.8) #encrusting HC
text(vct.scores[7,1],vct.scores[7,2]-0.02,
     row.names(vct.scores)[7],
     col='red', cex=0.8) #Foliose HC
text(vct.scores[8,1]-0.12,vct.scores[8,2]-0.05,
     row.names(vct.scores)[8],
     col='red', cex=0.8) #Massive HC
text(vct.scores[9,1]+0.09,vct.scores[9,2]-0.02,
     row.names(vct.scores)[9],
     col='red', cex=0.6) #Tabular HC
text(vct.scores[10,1]-0.1,vct.scores[10,2]+0.04,
     row.names(vct.scores)[10],
     col='red', cex=0.8) #Digitate OC
text(vct.scores[11,1],vct.scores[11,2]-0.02,
     row.names(vct.scores)[11],
     col='red', cex=0.8) #Lobate OC
text(vct.scores[12,1]+0.1,vct.scores[12,2]-0.05,
     row.names(vct.scores)[12],
     col='red', cex=0.6) #Massive OC
text(vct.scores[13,1],vct.scores[13,2]+0.02,
     row.names(vct.scores)[13],
     col='red', cex=0.6) #TWS
text(vct.scores[14,1]+0.03,vct.scores[14,2]+0.04,
     row.names(vct.scores)[14],
     col='red', cex=0.8)#turf
text(vct.scores[15,1]-0.1,vct.scores[15,2]-0.03,
     row.names(vct.scores)[15],
     col='red', cex=0.8) #Bushy OC
text(vct.scores[16,1],vct.scores[16,2],
     row.names(vct.scores)[16],
     col='red', cex=0.6) #can remove
text(vct.scores[18,1]-0.02,vct.scores[18,2]+0.02,
     row.names(vct.scores)[18],
     col='red', cex=0.8)#zonenc

text(vct.scores[20,1],vct.scores[20,2]+0.02,
     row.names(vct.scores)[20],
     col='red', cex=0.6) #Massive sponges
text(vct.scores[21,1]+0.02,vct.scores[21,2]+0.02,
     row.names(vct.scores)[21],
     col='red', cex=0.8) #massive zoanthids
text(vct.scores[22,1]-0.02,vct.scores[22,2]-0.1,
     row.names(vct.scores)[22],
     col='red', cex=0.55) #clustered oc
text(vct.scores[24,1]-0.05,vct.scores[24,2]+0.05,
     row.names(vct.scores)[24],
     col='red', cex=0.5)# others

text(vct.scores[26,1]-0.15,vct.scores[26,2]-0.01,
     row.names(vct.scores)[26],
     col='red', cex=0.6) #papillate sponge
text(vct.scores[28,1]+0.06,vct.scores[28,2]+0.13,
     row.names(vct.scores)[28],
     col='red', cex=0.8) #CCA
text(vct.scores[29,1],vct.scores[29,2]+0.03,
     row.names(vct.scores)[29],
     col='red', cex=0.8) #BSS
text(vct.scores[30,1]+0.15,vct.scores[30,2],
     row.names(vct.scores)[30],
     col='red', cex=0.8) #MA
text(vct.scores[31,1]+0.07,vct.scores[31,2],
     row.names(vct.scores)[31],
     col='red', cex=0.8) #US
groups2<-groups
groups2<-as.factor(groups2)
levels(groups2)<-my_cols
groups2<-as.vector(groups2)
points(comp_sub_RDA, display='site',
       scaling=3,
       col=groups2,
       pch=20, cex=1.5)
vct.scores1<-scores(comp_sub_RDA, 
                                choices=1:2, 
                                display='site',
                                scaling=3)
text(vct.scores1[3,1]+0.02,
     vct.scores1[3,2]-0.03,
     'Jialeshuei',
     col="#fc0f00",
     cex=0.6
     )
text(vct.scores1[20,1]-0.02,
     vct.scores1[20,2]+0.04,
     'Fenniaolin',
     col="#15e0fa",
     cex=0.6)
text(vct.scores1[22,1],
     vct.scores1[22,2]+0.03,
     'Bitou',
     col="#0c59fe" ,
     cex=0.6)
text(vct.scores1[21,1]+0.08,
     vct.scores1[21,2]+0.03,
     'Longdong',
     col="#0c59fe" ,
     cex=0.6)
legend('topleft', fill=my_cols,
       legend= c('East Coast','Kenting', 'Ludao', 'Lanyu', 'North Coast'), 
       title='Region', cex=1.2)

###PERMANOVA on benthic composition##################################
library(rempsyc)
comp_sub_std<-decostand(composition_substrates[,3:33], 'standardize')
comp_sub_hel_std<-decostand(comp_sub_hel, 'standardize')


comp_sub_hel_euc<-vegdist(comp_sub_hel_std, 'euc')


comp_sub_PERMANOVA2<-adonis2(comp_sub_hel_std~region, data=composition_substrates,
                            permutations = 999, method='euc')

comp_sub_PERMANOVA2
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = comp_sub_hel_std ~ region, data = composition_substrates, permutations = 999, method = "euc")
##          Df SumOfSqs     R2      F Pr(>F)    
## region    4   190.77 0.2564 1.7241  0.001 ***
## Residual 20   553.23 0.7436                  
## Total    24   744.00 1.0000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
comp_sub_PERMANOVA_nice<-nice_table(comp_sub_PERMANOVA2)
flextable::save_as_docx(comp_sub_PERMANOVA_nice, path = "Figures/Composition_PERMANOVA.docx")


#Pairwise Permutation
library(RVAideMemoire)
library(devtools)
#install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")
library(pairwiseAdonis)

permtest2<-pairwise.adonis(as.data.frame(comp_sub_hel_std), composition_substrates$region,
                           sim.method = 'euc')


permtest2<-as.data.frame(permtest2)
permtest2$sig<-NULL
permtest_nice2<-nice_table(permtest2,
                          note='* p < .05, ** p < .01, *** p < .001', highlight = T)
permtest_nice2

pairs

Df

SumsOfSqs

F.Model

R2

p

p.adjusted

KT vs LD

1.00

38.42

1.46

.15

.085

0.85

KT vs LY

1.00

35.30

1.28

.14

.099

0.99

KT vs EC

1.00

40.99

1.18

.13

.194

1.00

KT vs NC

1.00

66.34

2.18

.21

.009**

0.09

LD vs LY

1.00

32.68

1.59

.17

.054

0.54

LD vs EC

1.00

54.13

1.95

.20

.030*

0.30

LD vs NC

1.00

58.00

2.48

.24

.006**

0.06

LY vs EC

1.00

54.72

1.88

.19

.017*

0.17

LY vs NC

1.00

40.37

1.63

.17

.014*

0.14

EC vs NC

1.00

55.97

1.75

.18

.009**

0.09

Note. * p < .05, ** p < .01, *** p < .001

flextable::save_as_docx(permtest_nice2, path = "Figures/Composition_PairwisePermTest.docx")

6. Initial selection on benthic categories (sparsity + cross-validation)

##################CV by glmnet###########################
comp_sub_metrics<-comp_sub_metrics[,-26] #remove other motile invertebrates

#Set filter of sparsity = 0.5
library(glmnet)
## 載入需要的套件:Matrix
## 
## 載入套件:'Matrix'
## 下列物件被遮斷自 'package:tidyr':
## 
##     expand, pack, unpack
## Loaded glmnet 4.1-7
filter <- function(x, ...) which(colMeans(x == 0) > 0.5)
sparsity <- function(fraction = 0.7) {
  function(x, ...) which(colMeans(x == 0) > fraction)
}

#V4
cv_V4<-cv.glmnet(y=comp_sub_metrics$`VRM - 4 cm`, 
                 x=as.matrix(sqrt(comp_sub_metrics[,3:32])),
                 exclude=sparsity(0.5),
                 family=gaussian(),
                 nfolds=25)
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
V4_coef<-as.data.frame(as.matrix(coef(cv_V4, s='lambda.min')))

#PROC32
cv_PROC32<-cv.glmnet(y=comp_sub_metrics$`PROC - 32 cm`, 
                 x=as.matrix(sqrt(comp_sub_metrics[,3:32])),
                 exclude=sparsity(0.5),
                 family=Gamma('log'),
                 nfolds=25)
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
PROC32_coef<-as.data.frame(as.matrix(coef(cv_PROC32, s='lambda.min')))

#PLC4
cv_PLC4<-cv.glmnet(y=comp_sub_metrics$`PLC - 4 cm`, 
                   x=as.matrix(sqrt(comp_sub_metrics[,3:32])),
                   exclude=sparsity(0.5),
                   family=gaussian(),
                   nfolds=25)
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
PLC4_coef<-as.data.frame(as.matrix(coef(cv_PLC4, s='lambda.min')))

#PLC16
cv_PLC16<-cv.glmnet(y=comp_sub_metrics$`PLC - 16 cm`, 
                    x=as.matrix(sqrt(comp_sub_metrics[,3:32])),
                    exclude=sparsity(0.5),
                    family=Gamma('log'),
                    nfolds=25)
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
PLC16_coef<-as.data.frame(as.matrix(coef(cv_PLC16, s='lambda.min')))

#PLC32
cv_PLC32<-cv.glmnet(y=comp_sub_metrics$`PLC - 32 cm`, 
                    x=as.matrix(sqrt(comp_sub_metrics[,3:32])),
                    exclude=sparsity(0.5),
                    family=gaussian(),
                    nfolds=25)
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
PLC32_coef<-as.data.frame(as.matrix(coef(cv_PLC32, s='lambda.min')))

#D12
cv_D12<-cv.glmnet(y=comp_sub_metrics$`D (1 - 2 cm)`, 
                  x=as.matrix(sqrt(comp_sub_metrics[,3:32])),
                  exclude=sparsity(0.5),
                  family=Gamma('log'),
                  nfolds=25)
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
D12_coef<-as.data.frame(as.matrix(coef(cv_D12, s='lambda.min')))

#D24
cv_D24<-cv.glmnet(y=comp_sub_metrics$`D (2 - 4 cm)`, 
                  x=as.matrix(comp_sub_metrics[,3:32]),
                  exclude=sparsity(0.5),
                  family=Gamma('log'),
                  nfolds=25)
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
D24_coef<-as.data.frame(as.matrix(coef(cv_D24, s='lambda.min')))

#D816
cv_D816<-cv.glmnet(y=comp_sub_metrics$`D(8 - 16 cm)`, 
                   x=as.matrix(comp_sub_metrics[,3:32]),
                   exclude=sparsity(0.5),
                   family=Gamma('log'),
                   nfolds=25)
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
D816_coef<-as.data.frame(as.matrix(coef(cv_D816, s='lambda.min')))

#D1632
cv_D1632<-cv.glmnet(y=comp_sub_metrics$`D(16 - 32 cm)`, 
                    x=as.matrix(comp_sub_metrics[,3:32]),
                    exclude=sparsity(0.5),
                    family=gaussian(),
                    nfolds=25)
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
D1632_coef<-as.data.frame(as.matrix(coef(cv_D1632, s='lambda.min')))

#D3264
cv_D3264<-cv.glmnet(y=comp_sub_metrics$`D(32 - 64 cm)`, 
                    x=as.matrix(comp_sub_metrics[,3:32]),
                    exclude=sparsity(0.5),
                    family=gaussian(),
                    nfolds=25)
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
D3264_coef<-as.data.frame(as.matrix(coef(cv_D3264, s='lambda.min')))

#Sq
cv_VRSD<-cv.glmnet(y=comp_sub_metrics$Sq, 
                   x=as.matrix(comp_sub_metrics[,3:32]),
                   exclude=sparsity(0.5),
                   family=gaussian(),
                   nfolds=25)
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
Sq_coef<-as.data.frame(as.matrix(coef(cv_VRSD, s='lambda.min')))
#Show all plots
par(mfrow=c(2,3))
plot(cv_V4)
title('VRM - 4 cm', line= 3)
plot(cv_PLC4)
title('PLC - 4 cm', line = 3)
plot(cv_D12)
title('D (1 - 2 cm)', line = 3)
plot(cv_D24)
title('D (2 - 4 cm)', line = 3)
plot(cv_PLC16)
title('PLC - 16 cm', line = 3)
plot(cv_D816)
title('D (8 - 16 cm)', line = 3)

par(mfrow=c(2,3))
plot(cv_D1632)
title('D (16 - 32 cm)', line = 3)
plot(cv_PROC32)
title('PROC - 32 cm', line = 3)
plot(cv_PLC32)
title('PLC - 32 cm', line = 3)
plot(cv_D3264)
title('D (32 - 64 cm)', line = 3)
plot(cv_VRSD)
title('Sq', line = 3)

###Show all coefficients
cv_coef<-cbind(V4_coef, PLC4_coef, D12_coef, D24_coef,
      PLC16_coef, D816_coef, D1632_coef,
      PROC32_coef, PLC32_coef, D3264_coef, 
      Sq_coef)
colnames(cv_coef)<-c('VRM - 4 cm', 'PLC - 4 cm', 'D (1 - 2 cm)',
'D (2 - 4 cm)', 'PLC - 16 cm', 'D (8 - 16 cm)', 'D (16 - 32 cm)',
'PROC - 32 cm', 'PLC - 32 cm', 'D (32 - 64 cm)', 'Sq'     )
row.names(cv_coef)<-c('Intercept', 'Actinia', 'Ceratodictyon', 'Arborescent HC',
                      'Bushy HC', 'Columnar HC', 'Encrusting HC', 'Foliose HC', 
                      'Massive HC', 'Tabular HC', 'Digitate OC', 'Lobate OC', 
                      'Massive OC', 'TWS', 'Filamentous algae', 'Bushy OC', 
                      'Ascidian', 'Fan OC', 'Encrusting zoanthid', 'Globular sponge',
                      'Massive sponge', 'Massive zoanthid', 'Clustered OC', 'Encrusting OC',
                      'Repent sponge', 'Papillate sponge', 'Corallimorpharian', 'CCA', 
                      'BSS', 'Macroalgae', 'US')
cv_coef<-cv_coef[c(1,4,5,7,8,9,11,12,15,16,19,22,28, 29:31 ),]
categories_retained<-row.names(cv_coef)
cv_coef<-cbind(categories_retained,cv_coef)
colnames(cv_coef)[1]<-'Category'
cv_glm_nice<-nice_table(cv_coef)
cv_glm_nice

Category

VRM - 4 cm

PLC - 4 cm

D (1 - 2 cm)

D (2 - 4 cm)

PLC - 16 cm

D (8 - 16 cm)

D (16 - 32 cm)

PROC - 32 cm

PLC - 32 cm

D (32 - 64 cm)

Sq

Intercept

0.05

14.47

0.86

0.75

1.32

0.76

2.15

-0.94

1.98

2.23

0.33

Arborescent HC

0.05

0.00

0.08

0.36

0.00

0.00

0.00

0.81

0.00

0.00

0.00

Bushy HC

0.05

0.00

0.05

0.10

0.00

0.02

0.00

0.00

0.00

0.00

0.00

Encrusting HC

0.00

0.00

-0.11

-0.00

0.00

0.00

0.00

0.00

0.00

0.00

0.00

Foliose HC

0.00

0.00

-0.03

0.01

0.00

0.00

0.00

0.00

0.00

0.00

0.00

Massive HC

0.00

0.00

0.00

0.00

0.00

0.00

0.00

0.61

0.00

0.28

0.00

Digitate OC

0.00

0.00

-0.08

0.00

0.00

0.01

0.00

0.00

0.00

0.67

0.00

Lobate OC

0.05

0.00

0.04

0.00

0.00

0.00

0.00

0.43

0.00

0.00

0.00

Filamentous algae

0.03

0.00

-0.04

0.00

0.00

0.00

0.00

-0.57

0.00

-0.06

0.00

Bushy OC

0.00

0.00

0.02

0.44

0.00

0.00

0.00

0.00

0.00

0.00

0.00

Encrusting zoanthid

0.01

0.00

0.15

0.00

0.80

0.00

0.00

0.00

0.00

0.00

0.00

Massive zoanthid

0.00

0.00

-0.08

0.00

0.16

0.00

0.00

0.00

0.00

0.00

0.00

CCA

-0.03

0.00

-0.00

0.00

0.37

-0.13

-0.37

-1.44

0.00

-0.93

0.00

BSS

-0.00

0.00

-0.09

-0.00

0.00

0.00

-0.03

0.00

0.00

-0.02

0.00

Macroalgae

0.00

0.00

-0.04

0.00

0.00

0.00

0.00

0.07

0.00

0.00

0.00

US

0.00

0.00

-0.04

0.00

0.00

0.00

0.00

0.67

0.00

0.07

0.22

flextable::save_as_docx(cv_glm_nice, path = "Figures/CV_Coefficients.docx")

7. Second selection (stepwise VIF) and tb-RDA

Using stepwise VIF selection for the second selection on benthic categories.

Using hellinger-transformed composition to explain standardized complexity. Using permutation test to reveal significance. Using variance partitioning analysis to reveal the importance of each category.

par(mfrow=c(1,1))
library(dplyr)
#Standardize the hellinger-transformed composition
comp_sub_hel_std<-decostand(comp_sub_hel, 'standardize')

#Select the categories retained from initial selection
MetricsRDA_table<-dplyr::select(comp_sub_hel_std,`Arborescent HC`, `Bushy HC`, `Encrusting HC`,`Foliose HC`, `Massive HC`,
                         `Encrusting Zoanthids`,`Massive Zoanthids`, 
                         `Lobate OC`, `Digitate OC`,`Bushy OC`, `Filamentous algae`, `Macroalgae`, 
                         `CCA`,`BSS`, `US`
)

#Standardize the metrics
metrics_std<-decostand(metrics[,c(3,8,9,10,15,18,19,21,22,23,24)], 'standardize')
rda_table<- merge(MetricsRDA_table, metrics_std, by='row.names')
row.names(rda_table)<-rda_table$Row.names
rda_table<-rda_table[,-1]
rda_mod<-rda(rda_table[,c(16:26)]~., data=rda_table[,1:15])
rda_mod_summary<-summary(rda_mod)

#Second selection on benthic categories (stepwise VIF)
vif.cca(rda_mod)
##       `Arborescent HC`             `Bushy HC`        `Encrusting HC` 
##               2.759531               2.134464               6.744055 
##           `Foliose HC`           `Massive HC` `Encrusting Zoanthids` 
##               8.163998              14.994587              11.260313 
##    `Massive Zoanthids`            `Lobate OC`          `Digitate OC` 
##               5.858777               5.550216               4.022515 
##             `Bushy OC`    `Filamentous algae`             Macroalgae 
##               4.606060               3.227758               9.063365 
##                    CCA                    BSS                     US 
##               4.247677              19.328454               3.120272
#remove BSS 
#Select the retained categories again
MetricsRDA_table<-dplyr::select(comp_sub_hel_std,`Arborescent HC`, `Bushy HC`, `Encrusting HC`,`Foliose HC`, `Massive HC`,
                         `Encrusting Zoanthids`,`Massive Zoanthids`, 
                         `Lobate OC`, `Digitate OC`,`Bushy OC`, `Filamentous algae`, `Macroalgae`, 
                         `CCA`, `US`
)
metrics_std<-decostand(metrics[,c(3,8,9,10,15,18,19,21,22,23,24)], 'standardize')
rda_table<- merge(MetricsRDA_table, metrics_std, by='row.names')
row.names(rda_table)<-rda_table$Row.names
#Create a group, showing regions
groups<-c('NC', 'NC', 'LY', 'NC', 'LD', 'LD', 'LY', 'KT', 'EC', 'LD', 'LD',
          'KT', 'LY', 'EC', 'EC', 'KT', 'EC', 'EC', 'NC', 'KT', 'NC', 'LD',
          'LY', 'KT', 'LY')
groups2<-as.factor(groups)
levels(groups2)<-my_cols
rda_table<-rda_table[,-1]
rda_mod<-rda(rda_table[,c(15:25)]~., data=rda_table[,1:14])
rda_mod_summary<-summary(rda_mod)
vif.cca(rda_mod) #All VIF <10
##       `Arborescent HC`             `Bushy HC`        `Encrusting HC` 
##               2.599292               1.923347               5.332395 
##           `Foliose HC`           `Massive HC` `Encrusting Zoanthids` 
##               2.594163               3.685568               8.328774 
##    `Massive Zoanthids`            `Lobate OC`          `Digitate OC` 
##               5.541295               4.589677               3.215497 
##             `Bushy OC`    `Filamentous algae`             Macroalgae 
##               4.006500               3.174003               6.898325 
##                    CCA                     US 
##               2.960635               2.550567
##Visualization
#Show R^2
RsquareAdj(rda_mod)
## $r.squared
## [1] 0.7542331
## 
## $adj.r.squared
## [1] 0.4101595
par(mfrow=c(1,1))

ordiplot(rda_mod,type='n',
         scaling=3,,
         xlab=paste('RDA1, ', 
                    round(rda_mod_summary$concont$importance[2,1]*100, 2), '%'),
         ylab=paste('RDA2, ',round(rda_mod_summary$concont$importance[2,2]*100, 2), '%'),
         xlim=c(-1,1),
         cex=0.8)
title(
      sub='Variable selection by removing variables with VIF > 10')

rda_sc<-scores(rda_mod,choice=1:2,
               display='sp',
               scaling=2)


arrows(x0=0, y0=0,
       x1=rda_sc[,1],
       y1=rda_sc[,2],
       length=0.1,
       lty=1,
       lwd=1.5,
       col='red')
rda_sc1<-scores(rda_mod,choice=1:2,
               display='site',
               scaling=3)

points(x=rda_sc1[,1], y=rda_sc1[,2], pch=20, cex=1.5,
       col=as.vector(groups2))

ordihull(rda_mod, 
         groups= groups,
         draw='lines',
         scaling=3,
         alpha=0.5,
         lwd=1.5,
         col=my_cols,
         border = NULL)
rda_sc2<-scores(rda_mod,choice=1:2,
                scaling=2,
                display = 'bp')
row.names(rda_sc2)<-c('Arborescent HC', 'Bushy HC', 'Encrusting HC',
                      'Foliose HC', 'Massive HC',
                       'Encrusting Zoanthids','Massive Zoanthids', 
                       'Lobate OC', 'Digitate OC','Bushy OC', 'Filamentous algae', 
                      'Macroalgae','CCA', 'US')
arrows(x0=0, y0=0,
       x1=rda_sc2[,1],
       y1=rda_sc2[,2],
       length=0,
       lty=1,
       lwd=2,
       col='black')

text(x=rda_sc2[1,1]+0.15,
     y=rda_sc2[1,2]+0.065,
     row.names(rda_sc2)[1],
     cex=0.7,
     col='black')
text(x=rda_sc2[2,1]+0.265,
     y=rda_sc2[2,2]+0.1,
     row.names(rda_sc2)[2],
     cex=0.8,
     col='black')
text(x=rda_sc2[3,1]-0.13,
     y=rda_sc2[3,2]-0.12,
     row.names(rda_sc2)[3],
     cex=0.8,
     col='black')
text(x=rda_sc2[4,1]+0.25,
     y=rda_sc2[4,2]-0.07,
     row.names(rda_sc2)[4],
     cex=0.8,
     col='black')
text(x=rda_sc2[5,1]+0.18,
     y=rda_sc2[5,2]-0.2,
     row.names(rda_sc2)[5],
     cex=0.8,
     col='black')
text(x=rda_sc2[6,1]-0.46,
     y=rda_sc2[6,2]-0.017,
     row.names(rda_sc2)[6],
     cex=0.6,
     col='black')
text(x=rda_sc2[6,1]-0.3,
     y=rda_sc2[6,2]+0.05,
     row.names(rda_sc2)[7],
     cex=0.6,
     col='black')
text(x=rda_sc2[8,1]+0.37,
     y=rda_sc2[8,2]+0.07,
     row.names(rda_sc2)[8],
     cex=0.8,
     col='black')
text(x=rda_sc2[9,1]+0.35,
     y=rda_sc2[9,2],
     row.names(rda_sc2)[9],
     cex=0.8,
     col='black')
text(x=rda_sc2[10,1]+0.02,
     y=rda_sc2[10,2]+0.1,
     row.names(rda_sc2)[10],
     cex=0.8,
     col='black')
text(x=rda_sc2[11,1]-0.025,
     y=rda_sc2[11,2]+0.11,
     row.names(rda_sc2)[11],
     cex=0.6,
     col='black')
text(x=rda_sc2[12,1]+0.11,
     y=rda_sc2[12,2]-0.18,
     row.names(rda_sc2)[12],
     cex=0.55,
     col='black')
text(x=rda_sc2[13,1]-0.1,
     y=rda_sc2[13,2]+0.05,
     row.names(rda_sc2)[13],
     cex=0.8,
     col='black')
text(x=rda_sc2[14,1]+0.05,
     y=rda_sc2[14,2]-0.13,
     row.names(rda_sc2)[14],
     cex=0.8,
     col='black')
text(x=rda_sc[1,1]+0.3,
     y=rda_sc[1,2]+0.05,
     'VRM - 4 cm',
     cex=0.9,
     col='red')
text(x=rda_sc[2,1]-0.1,
     y=rda_sc[2,2]+0.25,
     'PLC - 4 cm',
     cex=0.9,
     col='red')
text(x=rda_sc[3,1]+0.45,
     y=rda_sc[3,2]-0.1,
     'PROC - 32 cm',
     cex=0.9,
     col='red')
text(x=rda_sc[4,1]+0.45,
     y=rda_sc[4,2]-0.2,
     'PLC - 32 cm',
     cex=0.8,
     col='red')
text(x=rda_sc[5,1]-0.4,
     y=rda_sc[5,2]-0.04,
     'PLC - 16 cm',
     cex=0.9,
     col='red')
text(x=rda_sc[6,1]+0.18,
     y=rda_sc[6,2]+0.2,
     'D (1 - 2 cm)',
     cex=0.9,
     col='red')
text(x=rda_sc[7,1]+0.15,
     y=rda_sc[7,2]+0.2,
     'D (2 - 4 cm)',
     cex=0.9,
     col='red')

text(x=rda_sc[8,1]+0.37,
     y=rda_sc[8,2]+0.05,
     'D (8 - 16 cm)',
     cex=0.9,
     col='red')
text(x=rda_sc[9,1]+0.5,
     y=rda_sc[9,2]-0.21,
     'D (16 - 32 cm)',
     cex=0.9,
     col='red')
text(x=rda_sc[10,1]+0.4,
     y=rda_sc[10,2]-0.31,
     'D (32 - 64 cm)',
     cex=0.9,
     col='red')
text(x=rda_sc[11,1]+0.15,
     y=rda_sc[11,2]-0.1,
     'Sq',
     cex=0.9,
     col='red')
legend('topright', fill=my_cols,
       legend= c('East Coast','Kenting', 'Ludao', 'Lanyu', 'North Coast'), 
       title='Region',
       cex=1.2)

#Show summary of the constrained model
rda_mod_summary 
## 
## Call:
## rda(formula = rda_table[, c(15:25)] ~ `Arborescent HC` + `Bushy HC` +      `Encrusting HC` + `Foliose HC` + `Massive HC` + `Encrusting Zoanthids` +      `Massive Zoanthids` + `Lobate OC` + `Digitate OC` + `Bushy OC` +      `Filamentous algae` + Macroalgae + CCA + US, data = rda_table[,      1:14]) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total          11.000     1.0000
## Constrained     8.297     0.7542
## Unconstrained   2.703     0.2458
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                         RDA1   RDA2   RDA3    RDA4    RDA5    RDA6     RDA7
## Eigenvalue            3.3298 1.8034 1.2039 1.00719 0.44074 0.22889 0.105438
## Proportion Explained  0.3027 0.1639 0.1094 0.09156 0.04007 0.02081 0.009585
## Cumulative Proportion 0.3027 0.4666 0.5761 0.66765 0.70772 0.72853 0.738116
##                           RDA8     RDA9    RDA10     RDA11    PC1    PC2
## Eigenvalue            0.086751 0.049569 0.034350 0.0066187 1.2629 0.6193
## Proportion Explained  0.007886 0.004506 0.003123 0.0006017 0.1148 0.0563
## Cumulative Proportion 0.746002 0.750509 0.753631 0.7542331 0.8690 0.9253
##                           PC3     PC4     PC5      PC6      PC7     PC8
## Eigenvalue            0.31843 0.21142 0.11673 0.091287 0.035766 0.02728
## Proportion Explained  0.02895 0.01922 0.01061 0.008299 0.003251 0.00248
## Cumulative Proportion 0.95429 0.97351 0.98412 0.992420 0.995672 0.99815
##                            PC9      PC10
## Eigenvalue            0.011693 0.0086399
## Proportion Explained  0.001063 0.0007854
## Cumulative Proportion 0.999215 1.0000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         RDA1   RDA2   RDA3   RDA4    RDA5    RDA6    RDA7
## Eigenvalue            3.3298 1.8034 1.2039 1.0072 0.44074 0.22889 0.10544
## Proportion Explained  0.4013 0.2174 0.1451 0.1214 0.05312 0.02759 0.01271
## Cumulative Proportion 0.4013 0.6187 0.7638 0.8852 0.93833 0.96592 0.97863
##                          RDA8     RDA9   RDA10     RDA11
## Eigenvalue            0.08675 0.049569 0.03435 0.0066187
## Proportion Explained  0.01046 0.005975 0.00414 0.0007978
## Cumulative Proportion 0.98909 0.995062 0.99920 1.0000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  4.03089 
## 
## 
## Species scores
## 
##           RDA1    RDA2     RDA3    RDA4      RDA5     RDA6
## V4      0.9858  0.3700  0.08551 -0.4348  0.048884 -0.13822
## PLC4   -0.4461  0.6915 -0.39389 -0.3200  0.136029  0.27246
## PROC32  0.8833 -0.4579 -0.06117  0.1519 -0.072617  0.12693
## PLC32   0.3694 -0.2888 -0.89377  0.3381  0.049175  0.08735
## PLC16  -0.3837 -0.1067 -0.25535 -0.4762 -0.688945  0.01808
## D1_2    0.7352  0.6647  0.23228  0.5269 -0.187704 -0.07290
## D2_4    0.6346  0.8562 -0.15508  0.1743 -0.224811  0.04522
## D8_16   0.8193  0.1599 -0.10335 -0.5615  0.204370  0.15831
## D16_32  0.7653 -0.4450 -0.29084 -0.3423  0.008021 -0.27649
## D32_64  0.6335 -0.5439  0.07099  0.2073 -0.128538  0.26853
## VRSD4   0.2863 -0.2234  0.75478 -0.2207 -0.060545  0.19929
## 
## 
## Site scores (weighted sums of species scores)
## 
##          RDA1     RDA2      RDA3     RDA4     RDA5      RDA6
## 82.5 -1.35323 -2.00650  1.930786 -0.43590 -1.82578  0.464597
## BT   -0.56245  0.17865 -0.287361 -1.74634 -1.55646  0.855593
## CCG   0.04190  1.73336  2.070078  0.04737  0.74665  0.602983
## CJ    0.07692 -0.73282  0.844272 -0.31376  1.71572  0.890174
## CK    0.63277 -1.32448  0.393205  0.71133 -0.42394 -0.006186
## DaBS -0.90423 -0.37032 -0.489845  0.20567  1.01853 -0.727456
## DC    0.23044  0.89952 -1.198009 -0.32731 -1.50122 -0.071492
## DiBS -0.77649  0.56325  0.574256  1.00468  0.93608 -0.735997
## FNL   1.08432  0.44154  0.005434 -0.12650  0.09923 -0.447854
## GGB   0.22928 -0.29771 -0.238655  0.07425 -0.33688  1.821578
## GW   -0.63263 -0.17139 -2.017465  0.89600  0.18060 -0.879072
## HBH   0.45026  0.03335  0.734766  0.89933 -0.50650  0.817287
## HR    0.73786  1.06863 -0.459736 -0.05137 -1.45310  0.278436
## HS    0.92370 -1.54142 -0.617109 -0.07919  0.45972  1.066678
## JH    0.36925 -0.70379  0.622378 -0.17409  1.32887  0.179194
## JLS   1.23257  0.68498  0.693771  0.33762 -0.03229  0.804830
## JMZ   1.59217 -0.70115 -1.046467 -0.84022  0.34913 -0.367680
## JQ    0.06906 -1.01478  0.607659  0.56581  2.01496  0.406128
## LD    0.05324  0.60557  1.113210 -2.57338  0.21174 -1.304019
## LDS  -0.08973  0.66957 -1.450940 -0.03823  0.73776  1.473230
## MYS  -2.17321  1.11339 -0.749820  0.15501  0.19930  0.688480
## SL    0.95862  0.17935  0.966558  1.71559 -0.81276 -0.918447
## TL   -1.04137  0.19287 -0.149560  0.32607 -0.02802 -2.339783
## TS    0.05251 -0.09990 -1.131960 -0.46018 -0.21007 -1.946109
## YY   -1.20152  0.60025 -0.719447  0.22775 -1.31129 -0.605095
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##          RDA1     RDA2      RDA3     RDA4     RDA5     RDA6
## 82.5 -0.84269 -2.04012  1.060994  0.01211 -1.04809  0.02771
## BT   -0.58393  0.23330 -0.247014 -1.79981 -1.66575  0.92760
## CCG   0.03141  1.89197  1.028257  0.40863  0.79207  0.30029
## CJ   -0.42795 -0.32892  0.402549 -0.78744  0.32720 -0.02688
## CK   -0.24026 -0.66382 -0.737456  0.76690 -0.96858 -0.62998
## DaBS -1.27160 -0.06655  0.399381  0.38701  0.99274 -0.16350
## DC    0.39397  0.81423 -0.588483  0.06418 -0.98809  0.32128
## DiBS -0.36121  0.08902  0.638646  0.57777  0.56759 -0.56014
## FNL   1.00879  0.55986 -0.003252 -0.10787 -0.03265 -0.61270
## GGB   0.29856 -0.30629 -0.056919 -0.29101 -0.39915  2.40790
## GW   -0.09683 -0.08946 -1.598913  0.84229  0.09702 -0.59836
## HBH   0.35403 -0.39126  0.645709  1.15982  0.19666  0.99702
## HR    0.59728  0.49384  0.459178  0.10309 -1.49311  0.05654
## HS    0.17828 -1.17547 -0.201842 -0.13487  0.44699 -0.25824
## JH    0.70265 -0.41194  1.345236 -0.34049  0.70598  0.61629
## JLS   1.02630  1.35732 -0.251769  0.12513 -0.21023  0.35993
## JMZ   1.71508 -0.74890 -0.804286 -0.67829  0.79514 -0.39024
## JQ   -0.20941 -0.71948 -0.241369  0.47293  1.10643  0.73041
## LD   -0.02623  0.43159  0.827087 -2.38006  0.80155 -1.29372
## LDS  -0.10623  0.39280 -1.828138  0.06004  0.69795  0.49406
## MYS  -1.76685  0.67139 -0.151129  0.17618  0.58388  0.60627
## SL    1.00704  0.05476  1.294795  1.49313 -0.68790 -0.97080
## TL   -0.80006  0.28761 -0.337091 -0.06179 -1.10194 -1.41786
## TS    0.69327 -0.87840 -0.829758 -0.36222  0.47264 -0.44032
## YY   -1.27339  0.54292 -0.224414  0.29465  0.01167 -0.48256
## 
## 
## Biplot scores for constraining variables
## 
##                            RDA1     RDA2       RDA3      RDA4      RDA5
## `Arborescent HC`        0.46608  0.32562 -0.1127738  0.184509 -0.398183
## `Bushy HC`              0.46969  0.42827 -0.1144653  0.218628  0.023679
## `Encrusting HC`        -0.06931 -0.25548 -0.0480181  0.051118 -0.351759
## `Foliose HC`            0.36195 -0.07125 -0.1436978  0.349544  0.033711
## `Massive HC`            0.32680 -0.28765  0.2218517 -0.067273 -0.471554
## `Encrusting Zoanthids` -0.19488  0.10841  0.1981194 -0.745849 -0.433235
## `Massive Zoanthids`    -0.19148  0.08691 -0.0003583 -0.621176 -0.416900
## `Lobate OC`             0.32639  0.10719 -0.0979112 -0.096798 -0.005932
## `Digitate OC`           0.20161  0.04163 -0.2664732 -0.226901  0.003099
## `Bushy OC`              0.27367  0.34742 -0.0588393  0.237920 -0.116191
## `Filamentous algae`    -0.20644  0.16162  0.3125838 -0.702297  0.206525
## Macroalgae              0.13039 -0.12046 -0.2565190 -0.177941  0.223476
## CCA                    -0.63980  0.39215  0.0171568 -0.145532 -0.514588
## US                      0.22531 -0.53728  0.2569424  0.004425  0.210593
##                             RDA6
## `Arborescent HC`        0.179640
## `Bushy HC`             -0.018169
## `Encrusting HC`        -0.147289
## `Foliose HC`            0.041287
## `Massive HC`           -0.076608
## `Encrusting Zoanthids` -0.037089
## `Massive Zoanthids`     0.148523
## `Lobate OC`            -0.002256
## `Digitate OC`           0.321167
## `Bushy OC`             -0.031960
## `Filamentous algae`    -0.177057
## Macroalgae              0.049085
## CCA                    -0.032177
## US                      0.290866
#Show the scores of RDA1 and RDA2
rda_scores<-as.data.frame(rda_sc2)
rda_scores<-round(rda_scores, 3)
rda_scores$Groups<-row.names(rda_scores)
rda_scores<-dplyr::select(rda_scores, Groups, RDA1, RDA2)
library(gt)
rda_scores_gt<-gt(rda_scores)
rda_scores_gt
Groups RDA1 RDA2
Arborescent HC 0.466 0.326
Bushy HC 0.470 0.428
Encrusting HC -0.069 -0.255
Foliose HC 0.362 -0.071
Massive HC 0.327 -0.288
Encrusting Zoanthids -0.195 0.108
Massive Zoanthids -0.191 0.087
Lobate OC 0.326 0.107
Digitate OC 0.202 0.042
Bushy OC 0.274 0.347
Filamentous algae -0.206 0.162
Macroalgae 0.130 -0.120
CCA -0.640 0.392
US 0.225 -0.537
gtsave(rda_scores_gt, 'Figures/RDA_scores.docx')
###############Check significance of RDA model##############
#Check the significance by axis
rda_aov<-anova.cca(rda_mod, by='axis',
                   step=5000)
rda_aov
## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = rda_table[, c(15:25)] ~ `Arborescent HC` + `Bushy HC` + `Encrusting HC` + `Foliose HC` + `Massive HC` + `Encrusting Zoanthids` + `Massive Zoanthids` + `Lobate OC` + `Digitate OC` + `Bushy OC` + `Filamentous algae` + Macroalgae + CCA + US, data = rda_table[, 1:14])
##          Df Variance       F Pr(>F)    
## RDA1      1   3.3298 16.0118  0.001 ***
## RDA2      1   1.8034  8.6720  0.110    
## RDA3      1   1.2039  5.7890  0.568    
## RDA4      1   1.0072  4.8433  0.628    
## RDA5      1   0.4407  2.1194  0.998    
## RDA6      1   0.2289  1.1007  1.000    
## RDA7      1   0.1054  0.5070  1.000    
## RDA8      1   0.0868  0.4172  1.000    
## RDA9      1   0.0496  0.2384  1.000    
## RDA10     1   0.0344  0.1652  1.000    
## RDA11     1   0.0066  0.0318  1.000    
## Residual 13   2.7034                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Significance at RDA1
rda_aov<-as.data.frame(rda_aov)
rda_aov$Axis<-row.names(rda_aov)
rda_aov<-dplyr::select(rda_aov, Axis, Df, Variance, `F`, `Pr(>F)`)
rda_aov<-nice_table(rda_aov, highlight = T)
rda_aov

Axis

Df

Variance

F

Pr(>F)

RDA1

1.00

3.33

16.01

0.00

RDA2

1.00

1.80

8.67

0.11

RDA3

1.00

1.20

5.79

0.57

RDA4

1.00

1.01

4.84

0.63

RDA5

1.00

0.44

2.12

1.00

RDA6

1.00

0.23

1.10

1.00

RDA7

1.00

0.11

0.51

1.00

RDA8

1.00

0.09

0.42

1.00

RDA9

1.00

0.05

0.24

1.00

RDA10

1.00

0.03

0.17

1.00

RDA11

1.00

0.01

0.03

1.00

Residual

13.00

2.70

flextable::save_as_docx(rda_aov, path = "Figures/RDA_ANOVA_ByAxis.docx")

#Check significance of the whole constrained model
rda_aov_all<-anova.cca(rda_mod,
                       step=5000) #Show overall significance
rda_aov_all
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = rda_table[, c(15:25)] ~ `Arborescent HC` + `Bushy HC` + `Encrusting HC` + `Foliose HC` + `Massive HC` + `Encrusting Zoanthids` + `Massive Zoanthids` + `Lobate OC` + `Digitate OC` + `Bushy OC` + `Filamentous algae` + Macroalgae + CCA + US, data = rda_table[, 1:14])
##          Df Variance      F Pr(>F)   
## Model    14   8.2966 2.1921  0.002 **
## Residual 10   2.7034                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rda_aov_all<-nice_table(rda_aov_all)
rda_aov_all

Df

Variance

F

Pr(>F)

14.00

8.30

2.19

0.00

10.00

2.70

flextable::save_as_docx(rda_aov_all, path = "Figures/RDA_ANOVA.docx")

#Significantly tells the variation
######################Variance partitioning###################################

varpart_plot<-varpart(rda_table[,c(15:25)],
                      X=rda_table[,1],
                      rda_table[,2],
                      rda_table[,3],
                      rda_table[,4])

varpart_table<-varpart_plot$part$fract[1:4,]
row.names(varpart_table)<-c('Arborescent HC','Bushy HC', 'Encrusting HC', 'Foliose HC' )
varpart_table<-round(varpart_table, 3)
varpart_plot<-varpart(rda_table[,c(15:25)],
                      X=rda_table[,5],
                      rda_table[,6],
                      rda_table[,7],
                      rda_table[,8])

varpart_table2<-varpart_plot$part$fract[1:4,]
row.names(varpart_table2)<-c('Massive HC','Encrusting zoanthids', 'Massive zoanthids', 'Lobate OC' )
varpart_table2<-round(varpart_table2, 3)
varpart_plot<-varpart(rda_table[,c(15:25)],
                      X=rda_table[,9],
                      rda_table[,10],
                      rda_table[,11],
                      rda_table[,12])

varpart_table3<-varpart_plot$part$fract[1:4,]
row.names(varpart_table3)<-c('Digitate OC','Bushy OC', 'Filamentous algae', 'Macroalgae' )
varpart_table3<-round(varpart_table3, 3)
varpart_plot<-varpart(rda_table[,c(15:25)],
                      rda_table[,13],
                      rda_table[,14] )

varpart_table4<-varpart_plot$part$fract[1:2,]
row.names(varpart_table4)<-c('CCA', 'US' )
varpart_table4<-round(varpart_table4, 3)
colnames(varpart_table4)[2:3]<-c('R.square', 'Adj.R.square')

varpart_all_group<-rbind(varpart_table, varpart_table2, varpart_table3, varpart_table4)
varpart_all_group$Testable<-NULL
varpart_all_group$Category<-row.names(varpart_all_group)
varpart_all_group<-dplyr::select(varpart_all_group, Category, Df, R.square, Adj.R.square)
library(gt)
varpart_gt<-gt(varpart_all_group)

varpart_gt
Category Df R.square Adj.R.square
Arborescent HC 1 0.096 0.056
Bushy HC 1 0.104 0.065
Encrusting HC 1 0.023 -0.020
Foliose HC 1 0.058 0.017
Massive HC 1 0.062 0.021
Encrusting zoanthids 1 0.077 0.036
Massive zoanthids 1 0.056 0.015
Lobate OC 1 0.039 -0.003
Digitate OC 1 0.030 -0.012
Bushy OC 1 0.052 0.011
Filamentous algae 1 0.076 0.036
Macroalgae 1 0.022 -0.020
CCA 1 0.162 0.125
US 1 0.077 0.037
gtsave(varpart_gt, 'Figures/Varpart.docx')

8. Displaying scatter plot to examine the relationships between complexity metrics and each category

library(glmmTMB)
library(buildmer)
library(lme4)
## 
## 載入套件:'lme4'
## 下列物件被遮斷自 'package:RVAideMemoire':
## 
##     dummy
## 下列物件被遮斷自 'package:raster':
## 
##     getData
library(ggplot2)
library(tidyr)

#Select categories with sparsity < 0.5


f<-which(colSums(comp_sub_metrics[,3:32]==0)<0.5*nrow(comp_sub_metrics)) #Select categories with sparsity < 0.5

colnames(comp_sub_metrics)[33:56]<-c('S4', 'T4', 'V4', 'S32', 'T32',
                                     'V32', 'PROC4', 'PLC4', 'PROC32',
                                     'PLC32', 'S16', 'T16', 'V16', 'PROC16',
                                     'PLC16', 'D64', 'SC', 'D1_2', 'D2_4', 'D4_8', 'D8_16',
                                     'D16_32', 'D32_64', 'Sq')

glmm_table1<-comp_sub_metrics[,c(2,f+2, 35)]
glmm_table1$region<-as.factor(glmm_table1$region)
col_names<-c('Region', 'Arborescent HC','Bushy HC', 'Encrusting HC',
  'Foliose HC', 'Massive HC', 'Digitate OC',
  'Lobate OC', 'Filamentous algae', 'Bushy OC',
  'Encrusting zoanthids', 'Massive zoanthids', 'CCA', 'BSS', 'Macroalgae',
  'US')



#V4
table<-pivot_longer(data=glmm_table1, cols=-c('V4','region'), names_to = 'MFG', values_to ='cover' )
V4_scater<-ggplot(data=table, aes(y=V4,x=cover))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  facet_wrap(vars(MFG), scales='free')+
  scale_color_manual(values=my_cols) + 
  theme(axis.text=element_text(size=10))


#PROC32
glmm_table2<-comp_sub_metrics[,c(2,f+2, 41)]
glmm_table2[2:15]<-glmm_table2[2:15]
#colnames(glmm_table2)[1:16]<-col_names

table<-pivot_longer(data=glmm_table2, cols=-c('PROC32','region'), names_to = 'MFG', values_to ='cover' )
PROC32_scater<-ggplot(data=table, aes(y=PROC32, x=cover))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  facet_wrap(vars(MFG), scales='free')+
  scale_color_manual(values=my_cols)+ 
  theme(axis.text=element_text(size=10))

#PLC4
glmm_table12<-comp_sub_metrics[,c(2,f+2, 40)]
glmm_table12[2:15]<-glmm_table12[2:15]

table<-pivot_longer(data=glmm_table12, cols=-c('PLC4','region'), names_to = 'MFG', values_to ='cover' )
PLC4_scater<-ggplot(data=table, aes(y=PLC4, x=cover))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  facet_wrap(vars(MFG), scales='free')+
  scale_color_manual(values=my_cols)+ 
  theme(axis.text=element_text(size=10))

#PLC16
glmm_table4<-comp_sub_metrics[,c(2,f+2, 50-2-1)]
glmm_table4[2:15]<-glmm_table4[2:15]
table<-pivot_longer(data=glmm_table4, cols=-c('PLC16','region'), names_to = 'MFG', values_to ='cover' )
PLC16_scater<-ggplot(data=table, aes(y=PLC16, x=cover))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  facet_wrap(vars(MFG), scales='free')+
  scale_color_manual(values=my_cols)+ 
  theme(axis.text=element_text(size=10))

#PLC32
glmm_table3<-comp_sub_metrics[,c(2,f+2, 45-2-1)]
glmm_table3[2:15]<-glmm_table3[2:15]
table<-pivot_longer(data=glmm_table3, cols=-c('PLC32','region'), names_to = 'MFG', values_to ='cover' )
PLC32_scater<-ggplot(data=table, aes(y=PLC32, x=cover))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  facet_wrap(vars(MFG), scales='free')+
  scale_color_manual(values=my_cols)+ 
  theme(axis.text=element_text(size=10))

#D12
glmm_table5<-comp_sub_metrics[,c(2,f+2, 53-2-1)]
glmm_table5[2:15]<-glmm_table5[2:15]
table<-pivot_longer(data=glmm_table5, cols=-c('D1_2','region'), names_to = 'MFG', values_to ='cover' )
D12_scater<-ggplot(data=table, aes(y=D1_2, x=cover))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  facet_wrap(vars(MFG), scales='free')+
  scale_color_manual(values=my_cols)+ 
  theme(axis.text=element_text(size=10))

#D2_4
glmm_table6<-comp_sub_metrics[,c(2,f+2, 54-2-1)]
glmm_table6[2:15]<-glmm_table6[2:15]
table<-pivot_longer(data=glmm_table6, cols=-c('D2_4','region'), names_to = 'MFG', values_to ='cover' )
D24_scater<-ggplot(data=table, aes(y=`D2_4`, x=cover ))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  facet_wrap(vars(MFG), scales='free')+
  scale_color_manual(values=my_cols)+ 
  theme(axis.text=element_text(size=10))

#D8_16
glmm_table8<-comp_sub_metrics[,c(2,f+2, 56-2-1)]
glmm_table8[2:15]<-glmm_table8[2:15]
table<-pivot_longer(data=glmm_table8, cols=-c('D8_16','region'), names_to = 'MFG', values_to ='cover' )
D816_scater<-ggplot(data=table, aes(y=`D8_16`, x=cover))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  facet_wrap(vars(MFG), scales='free')+
  scale_color_manual(values=my_cols)+ 
  theme(axis.text=element_text(size=10))

#D16_32
glmm_table9<-comp_sub_metrics[,c(2,f+2, 57-2-1)]
glmm_table9[2:15]<-glmm_table9[2:15]
table<-pivot_longer(data=glmm_table9, cols=-c('D16_32','region'), names_to = 'MFG', values_to ='cover' )
D1632_scater<-ggplot(data=table, aes(y=`D16_32`, x=cover))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  facet_wrap(vars(MFG), scales='free')+
  scale_color_manual(values=my_cols)+ 
  theme(axis.text=element_text(size=10))

#D32_64
glmm_table10<-comp_sub_metrics[,c(2,f+2, 58-2-1)]
glmm_table10[2:15]<-glmm_table10[2:15]
table<-pivot_longer(data=glmm_table10, cols=-c('D32_64','region'), names_to = 'MFG', values_to ='cover' )
D3264_scater<-ggplot(data=table, aes(y=`D32_64`, x=cover))+
  geom_point(aes(col=region))+
  geom_smooth(method='lm')+
  facet_wrap(vars(MFG), scales='free')+
  scale_color_manual(values=my_cols)+ 
  theme(axis.text=element_text(size=10))

#Sq
glmm_table11<-comp_sub_metrics[,c(2,f+2, 59-2-1)]
glmm_table11[2:15]<-glmm_table11[2:15]
table<-pivot_longer(data=glmm_table11, cols=-c('Sq','region'), names_to = 'MFG', values_to ='cover' )
Sq_scater<-ggplot(data=table, aes(y=Sq, x=cover))+
  geom_point(aes( col=region))+
  geom_smooth(method='lm')+
  facet_wrap(vars(MFG), scales='free')+
  scale_color_manual(values=my_cols)+ 
  theme(axis.text=element_text(size=10))


#Show all plots
V4_scater
## `geom_smooth()` using formula = 'y ~ x'

PROC32_scater
## `geom_smooth()` using formula = 'y ~ x'

PLC4_scater
## `geom_smooth()` using formula = 'y ~ x'

PLC16_scater
## `geom_smooth()` using formula = 'y ~ x'

PLC32_scater
## `geom_smooth()` using formula = 'y ~ x'

D12_scater
## `geom_smooth()` using formula = 'y ~ x'

D24_scater
## `geom_smooth()` using formula = 'y ~ x'

D816_scater
## `geom_smooth()` using formula = 'y ~ x'

D1632_scater
## `geom_smooth()` using formula = 'y ~ x'

D3264_scater
## `geom_smooth()` using formula = 'y ~ x'

Sq_scater
## `geom_smooth()` using formula = 'y ~ x'

9.1. Confirmation of relationships by using GLMM

The coefficients were visualized by dwplots.

########################################GLMM########################################
#################################vif prior to modelling######################################
library(dplyr)
library(usdm)
glmm_vif_table<-dplyr::select(glmm_table11,hcarb, hcbus,hcenc,hcfol,
                       hcmas, zoaenc,zoamas, ocdig,oclob, ocbus,ma, fil, cru_or_spoenc, us)
vifstep(glmm_vif_table) #oclob and zoaenc are removed


###################GLMM    modelling###################################################
#####################################################################
colnames(glmm_table1)[17]<-'V4'
V4_glmm<-buildglmmTMB(V4~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
                        fil+cru_or_spoenc+ma+us+(1|region), 
                      data=glmm_table1, family=gaussian(),
                      buildmerControl = list(include = ~ (1|region), direction='forward'))

V4_mod_sum<-summary(V4_glmm)



colnames(glmm_table12)[c(1,17)]<-c('region','PLC4')

PLC4_glmm<-buildglmmTMB(PLC4~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
                          fil+cru_or_spoenc+ma+us+(1|region), 
                        data=glmm_table12, family=gaussian(),
                        buildmerControl = list(include = ~ (1|region), direction='forward'))
PLC4_mod_sum<-summary(PLC4_glmm)
#plot(PLC4_glmm)


colnames(glmm_table4)[17]<-'PLC16'
PLC16_glmm<-buildglmmTMB(PLC16~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
                           fil+cru_or_spoenc+ma+us+(1|region), 
                         data=glmm_table4, family=Gamma('log'),
                         buildmerControl = list(include = ~ (1|region), direction='forward'))
PLC16_mod_sum<-summary(PLC16_glmm)
#plot(PLC32_glmm)


colnames(glmm_table3)[17]<-'PLC32'
PLC32_glmm<-buildglmmTMB(PLC32~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
                           fil+cru_or_spoenc+ma+us+(1|region), 
                         data=glmm_table3, family=gaussian(),
                         buildmerControl = list(include = ~ (1|region), direction='forward'))
PLC32_mod_sum<-summary(PLC32_glmm)
#plot(PLC32_glmm)


colnames(glmm_table5)[17]<-'D1_2'
D12_glmm<-buildglmmTMB(D1_2~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
                         fil+cru_or_spoenc+ma+us+(1|region), 
                       data=glmm_table5, family=Gamma('log'),
                       buildmerControl = list(include = ~ (1|region), direction='forward'))
D12_mod_sum<-summary(D12_glmm)
#plot(D12_glmm)


colnames(glmm_table6)[17]<-'D2_4'
D24_glmm<-buildglmmTMB(D2_4~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
                         fil+cru_or_spoenc+ma+us+(1|region), 
                       data=glmm_table6, family=Gamma('log'),
                       buildmerControl = list(include = ~ (1|region), direction='forward'))
D24_mod_sum<-summary(D24_glmm)




colnames(glmm_table8)[17]<-'D8_16'
D816_glmm<-buildglmmTMB(D8_16~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
                          fil+cru_or_spoenc+ma+us+(1|region), 
                        data=glmm_table8, family=Gamma('log'),
                        buildmerControl = list(include = ~ (1|region), direction='forward'))
D816_mod_sum<-summary(D816_glmm)


colnames(glmm_table9)[17]<-'D16_32'
D1632_glmm<-buildglmmTMB(D16_32~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
                           fil+cru_or_spoenc+ma+us+(1|region), 
                         data=glmm_table9, family=gaussian(),
                         buildmerControl = list(include = ~ (1|region), direction='forward'))
D1632_mod_sum<-summary(D1632_glmm)
#plot(D1632_glmm)


colnames(glmm_table10)[17]<-'D32_64'
D3264_glmm<-buildglmmTMB(D32_64~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
                           fil+cru_or_spoenc+ma+us+(1|region), 
                         data=glmm_table10, family=gaussian(),
                         buildmerControl = list(include = ~ (1|region), direction='forward'))
D3264_mod_sum<-summary(D3264_glmm)


colnames(glmm_table11)[17]<-'Sq'
VRSD_glmm<-buildglmmTMB(Sq~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
                          fil+cru_or_spoenc+ma+us+(1|region), 
                        data=glmm_table11, family=gaussian(),
                        buildmerControl = list(include = ~ (1|region), direction='forward'))
VRSD_mod_sum<-summary(VRSD_glmm) # No benthic variables can significantly explain Sq

colnames(glmm_table2)[17]<-'PROC32'
PROC32_glmm<-buildglmmTMB(PROC32~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
                            fil+cru_or_spoenc+ma+us+(1|region), 
                          data=glmm_table2, family=Gamma('log'),
                          buildmerControl = list(include = ~ (1|region)))

PROC32_mod_sum<-summary(PROC32_glmm)
###########Show coefficients#########################

library(broom.mixed)
library(dotwhisker)

V4_mod<-glmmTMB(V4_mod_sum$call$formula,data=glmm_table1,family=gaussian())
V4_mod_sum<-summary(V4_mod)

V4_coef<-dwplot(V4_mod, effects='fixed',whisker_args = list(color = "blue"),
                dot_args = list(color = "black") )+
  geom_vline(xintercept=0,lty=2)+
  scale_y_discrete(labels=c('hcbus'='Bushy HC',
                            'fil'='Filamentous algae',
                            'ocbus' = 'Bushy OC',
                            'cru_or_spoenc' = 'CCA',
                            'hcarb' = 'Arborescent HC',
                            'zoamas' = 'Massive zoanthids'))+
  ggtitle('VRM - 4 cm')

V4_coef

PROC32_mod<-glmmTMB(PROC32_mod_sum$call$formula,data=glmm_table2,family=Gamma('log'))
PROC32_mod_sum<-summary(PROC32_mod)
PROC32_coef<-dwplot(PROC32_mod, effects='fixed',whisker_args = list(color = "red"),
                    dot_args = list(color = "black"))+geom_vline(xintercept=0,lty=2)+
 scale_y_discrete(labels=c('cru_or_spoenc' = 'CCA',
                           'hcmas' = 'Massive HC',
                           'fil' = 'Filamentous algae'))+
  ggtitle('PROC - 32 cm')
PROC32_coef

PLC4_mod<-glmmTMB(PLC4_mod_sum$call$formula,data=glmm_table12,family=gaussian())

PLC4_mod_sum<-summary(PLC4_mod)

PLC4_coef<-dwplot(PLC4_mod, effects='fixed',whisker_args = list(color = "blue"),
                  dot_args = list(color = "black") )+geom_vline(xintercept=0,lty=2)+
  scale_y_discrete(labels= c('hcmas'='Massive HC','cru_or_spoenc'='CCA'))+
  ggtitle('PLC - 4 cm')
PLC4_coef

PLC16_mod<-glmmTMB(PLC16_mod_sum$call$formula,data=glmm_table4,family=Gamma('log'))

PLC16_mod_sum<-summary(PLC16_mod)
PLC16_coef<-dwplot(PLC16_mod, effects='fixed',whisker_args = list(color = "orange"),
                   dot_args = list(color = "black"))+geom_vline(xintercept=0,lty=2)+
scale_y_discrete(labels=c('cru_or_spoenc' = 'CCA',
                          'ma' = 'Macroalgae'))+
  ggtitle('PLC - 16 cm')
PLC16_coef

PLC32_mod<-glmmTMB(PLC32_mod_sum$call$formula,data=glmm_table3,family=gaussian())
PLC32_mod_sum<-summary(PLC32_mod)
PLC32_coef<-dwplot(PLC32_mod, effects='fixed',whisker_args = list(color = "red"),
                   dot_args = list(color = "black"))+geom_vline(xintercept=0,lty=2)+
  scale_y_discrete(labels=c('fil' = 'Filamentous algae'))+
  ggtitle('PLC - 32 cm')
PLC32_coef

D12_mod<-glmmTMB(D12_mod_sum$call$formula,data=glmm_table5,family=Gamma('log'))
D12_mod_sum<-summary(D12_mod)
D12_coef<-dwplot(D12_mod, effects='fixed',whisker_args = list(color = "blue"),
                 dot_args = list(color = "black"))+geom_vline(xintercept=0,lty=2)+
  scale_y_discrete(labels=c('hcarb' = 'Arborescent HC',
                            'hcbus' = 'Bushy HC',
                            'ocbus' = 'Bushy OC',
                            'hcfol' = 'Foliose HC',
                            'hcenc' = 'Encrusting HC',
                            'us' = 'US',
                            'zoamas' = 'Massive zoanthids'))+
  ggtitle('D (1 - 2 cm)')
D12_coef

D24_mod<-glmmTMB(D24_mod_sum$call$formula,data=glmm_table6,family=Gamma('log'))
D24_mod_sum<-summary(D24_mod)
D24_coef<-dwplot(D24_mod, effects='fixed',whisker_args = list(color = "blue"),
       dot_args = list(color = "black"))+geom_vline(xintercept=0,lty=2)+
  scale_y_discrete(labels=c('hcarb' = 'Arborescent HC',
                            'hcbus' = 'Bushy HC',
                            'ocbus' = 'Bushy OC',
                            'hcfol' = 'Foliose HC',
                            'fil' = 'Filamentous algae',
                            'cru_or_spoenc' = 'CCA'))+
  ggtitle('D (2 - 4 cm)')
D24_coef

D816_mod<-glmmTMB(D816_mod_sum$call$formula,data=glmm_table8,family=Gamma('log'))
D816_mod_sum<-summary(D816_mod)
D816_coef<-dwplot(D816_mod, effects='fixed',whisker_args = list(color = "orange"),
       dot_args = list(color = "black"))+geom_vline(xintercept=0,lty=2)+
  scale_y_discrete(labels=c('cru_or_spoenc' = 'CCA',
                            'zoamas' = 'Massive zoanthids',
                            'hcfol' = 'Foliose HC',
                            'fil' = 'Filamentous algae',
                            'ocdig' = 'Digitate OC'))+
  ggtitle('D (8 - 16 cm)')
D816_coef

D1632_mod<-glmmTMB(D1632_mod_sum$call$formula,data=glmm_table9,family=gaussian())
D1632_mod_sum<-summary(D1632_mod)
D1632_coef<-dwplot(D1632_mod, effects='fixed',whisker_args = list(color = "orange"),
       dot_args = list(color = "black"))+geom_vline(xintercept=0,lty=2)+
 scale_y_discrete(labels=c('cru_or_spoenc' = 'CCA'))+
  ggtitle('D (16 - 32 cm)')
D1632_coef

D3264_mod<-glmmTMB(D3264_mod_sum$call$formula,data=glmm_table10,family=gaussian())
D3264_mod_sum<-summary(D3264_mod)
D3264_coef<-dwplot(D3264_mod, effects='fixed',whisker_args = list(color = "red"),
                   dot_args = list(color = "black"))+geom_vline(xintercept=0,lty=2)+
 scale_y_discrete(labels=c('cru_or_spoenc' = 'CCA',
                           'hcmas' = 'Massive HC'))+
  ggtitle('D (32 - 64 cm)')
D3264_coef

VRSD_mod<-glmmTMB(VRSD_mod_sum$call$formula,data=glmm_table11,family=gaussian())

9.2. Showing the mean value and 95% CI of each coefficient in GLMMs

####################Show CI and coefficients of every model#################################
#V4
V4_CI<-confint(V4_mod)[-nrow(confint(V4_mod)),]
stats_table <- as.data.frame(V4_mod_sum$coefficients$cond)
stats_table_V4<-cbind(
  row.names(stats_table),
  round(stats_table, 3), V4_CI[,-3]
)
colnames(stats_table_V4) <- c(
  "Term", "B", "SE", "t", "p",
  "CI_lower", "CI_upper")
stats_table_V4$Term<-c('Intercept','Bushy HC', 'Filamentous algae', 
        'Bushy OC', 'CCA', 'Arborescent HC', 'Massive Zoanthids' )
V4_nt<-nice_table(stats_table_V4, highlight = TRUE,
                  note="* p < .05, ** p < .01, *** p < .001")
V4_nt

Term

β

SE

t

p

95% CI

Intercept

0.050

0.00

12.36

< .001***

[0.04, 0.06]

Bushy HC

0.193

0.04

5.33

< .001***

[0.12, 0.26]

Filamentous algae

0.103

0.02

4.98

< .001***

[0.06, 0.14]

Bushy OC

0.673

0.19

3.60

< .001***

[0.31, 1.04]

CCA

-0.219

0.07

-3.20

.001**

[-0.35, -0.09]

Arborescent HC

0.345

0.15

2.36

.018*

[0.06, 0.63]

Massive Zoanthids

0.075

0.03

2.22

.026*

[0.01, 0.14]

Note. * p < .05, ** p < .01, *** p < .001

#PROC32
PROC32_CI<-confint(PROC32_mod)[-nrow(confint(PROC32_mod)),]
stats_table <- as.data.frame(PROC32_mod_sum$coefficients$cond)
stats_table_PROC32<-cbind(
  row.names(stats_table),
  round(stats_table,3), PROC32_CI[,-3]
)
colnames(stats_table_PROC32) <- c(
  "Term", "B", "SE", "t", "p",
  "CI_lower", "CI_upper")
stats_table_PROC32$Term<-c('Intercept','CCA', 'Massive HC', 
                       'Filamentous algae' )
PROC32_nt<-nice_table(stats_table_PROC32, highlight = TRUE,
                      note="* p < .05, ** p < .01, *** p < .001")
PROC32_nt

Term

β

SE

t

p

95% CI

Intercept

-0.755

0.10

-7.40

< .001***

[-0.96, -0.56]

CCA

-7.450

1.96

-3.80

< .001***

[-11.30, -3.60]

Massive HC

2.008

0.81

2.48

.013*

[0.42, 3.59]

Filamentous algae

-1.411

0.50

-2.80

.005**

[-2.40, -0.42]

Note. * p < .05, ** p < .01, *** p < .001

#PLC4
PLC4_CI<-confint(PLC4_mod)[-nrow(confint(PLC4_mod)),]
stats_table <- as.data.frame(PLC4_mod_sum$coefficients$cond)
stats_table_PLC4<-cbind(
  row.names(stats_table),
  round(stats_table,3), PLC4_CI[,-3]
)
colnames(stats_table_PLC4) <- c(
  "Term", "B", "SE", "t", "p",
  "CI_lower", "CI_upper")
stats_table_PLC4$Term<-c('Intercept','Massive HC', 
                           'CCA' )
PLC4_nt<-nice_table(stats_table_PLC4, highlight = TRUE,
                    note="* p < .05, ** p < .01, *** p < .001")
PLC4_nt

Term

β

SE

t

p

95% CI

Intercept

14.469

0.75

19.39

< .001***

[13.01, 15.93]

Massive HC

-19.127

6.40

-2.99

.003**

[-31.68, -6.57]

CCA

46.577

15.81

2.95

.003**

[15.60, 77.56]

Note. * p < .05, ** p < .01, *** p < .001

#PLC16
PLC16_CI<-confint(PLC16_mod)[-nrow(confint(PLC16_mod)),]
stats_table <- as.data.frame(PLC16_mod_sum$coefficients$cond)
stats_table_PLC16<-cbind(
  row.names(stats_table),
  round(stats_table,3), PLC16_CI[,-3]
)
colnames(stats_table_PLC16) <- c(
  "Term", "B", "SE", "t", "p",
  "CI_lower", "CI_upper")
stats_table_PLC16$Term<-c('Intercept','CCA', 
                         'Macroalgae' )
PLC16_nt<-nice_table(stats_table_PLC16, highlight = TRUE,
                     note="* p < .05, ** p < .01, *** p < .001")
PLC16_nt

Term

β

SE

t

p

95% CI

Intercept

1.305

0.05

26.62

< .001***

[1.21, 1.40]

CCA

3.324

1.11

3.00

.003**

[1.15, 5.50]

Macroalgae

0.277

0.48

0.58

.565

[-0.66, 1.22]

Note. * p < .05, ** p < .01, *** p < .001

#PLC32
PLC32_CI<-confint(PLC32_mod)[-nrow(confint(PLC32_mod)),]
stats_table <- as.data.frame(PLC32_mod_sum$coefficients$cond)
stats_table_PLC32<-cbind(
  row.names(stats_table),
  round(stats_table,3), PLC32_CI[,-3]
)
colnames(stats_table_PLC32) <- c(
  "Term", "B", "SE", "t", "p",
  "CI_lower", "CI_upper")
stats_table_PLC32$Term<-c('Intercept','Filamentous algae' )
PLC32_nt<-nice_table(stats_table_PLC32,   highlight=T, note="* p < .05, ** p < .01, *** p < .001")
PLC32_nt

Term

β

SE

t

p

95% CI

Intercept

2.056

0.09

22.68

< .001***

[1.88, 2.23]

Filamentous algae

-2.633

0.94

-2.79

.005**

[-4.48, -0.78]

Note. * p < .05, ** p < .01, *** p < .001

#D12
D12_CI<-confint(D12_mod)[-nrow(confint(D12_mod)),]
stats_table <- as.data.frame(D12_mod_sum$coefficients$cond)
stats_table_D12<-cbind(
  row.names(stats_table),
  round(stats_table,3), D12_CI[,-3]
)
colnames(stats_table_D12) <- c(
  "Term", "B", "SE", "t", "p",
  "CI_lower", "CI_upper")
stats_table_D12$Term<-c('Intercept','Arborescent HC','Bushy HC', 'Bushy OC',
                        'Foliose HC', 'Encrusting HC', 'US', 'Massive zoanthids')
D12_nt<-nice_table(stats_table_D12,   highlight=T, note="* p < .05, ** p < .01, *** p < .001")
D12_nt

Term

β

SE

t

p

95% CI

Intercept

0.753

0.01

108.26

< .001***

[0.74, 0.77]

Arborescent HC

0.566

0.19

2.98

.003**

[0.19, 0.94]

Bushy HC

0.134

0.04

3.03

.002**

[0.05, 0.22]

Bushy OC

0.484

0.23

2.07

.039*

[0.02, 0.94]

Foliose HC

0.058

0.03

1.89

.059

[-0.00, 0.12]

Encrusting HC

-0.039

0.02

-1.75

.080

[-0.08, 0.00]

US

-0.053

0.05

-1.09

.274

[-0.15, 0.04]

Massive zoanthids

-0.021

0.04

-0.50

.616

[-0.10, 0.06]

Note. * p < .05, ** p < .01, *** p < .001

#D24
D24_CI<-confint(D24_mod)[-nrow(confint(D24_mod)),]
stats_table <- as.data.frame(D24_mod_sum$coefficients$cond)
stats_table_D24<-cbind(
  row.names(stats_table),
  round(stats_table,3), D24_CI[,-3]
)
colnames(stats_table_D24) <- c(
  "Term", "B", "SE", "t", "p",
  "CI_lower", "CI_upper")
stats_table_D24$Term<-c('Intercept','Arborescent HC','Bushy HC', 'Bushy OC',
                        'Foliose HC', 'Filamentous algae', 'CCA')
D24_nt<-nice_table(stats_table_D24,   highlight=T, note="* p < .05, ** p < .01, *** p < .001")
D24_nt

Term

β

SE

t

p

95% CI

Intercept

0.738

0.00

223.45

< .001***

[0.73, 0.74]

Arborescent HC

0.520

0.11

4.63

< .001***

[0.30, 0.74]

Bushy HC

0.152

0.03

5.58

< .001***

[0.10, 0.20]

Bushy OC

0.721

0.14

5.09

< .001***

[0.44, 1.00]

Foliose HC

0.051

0.01

3.27

.001**

[0.02, 0.08]

Filamentous algae

0.041

0.02

2.59

.010*

[0.01, 0.07]

CCA

0.110

0.05

2.12

.034*

[0.01, 0.21]

Note. * p < .05, ** p < .01, *** p < .001

#D816
D816_CI<-confint(D816_mod)[-nrow(confint(D816_mod)),]
stats_table <- as.data.frame(D816_mod_sum$coefficients$cond)
stats_table_D816<-cbind(
  row.names(stats_table),
  round(stats_table,3), D816_CI[,-3]
)
colnames(stats_table_D816) <- c(
  "Term", "B", "SE", "t", "p",
  "CI_lower", "CI_upper")
stats_table_D816$Term<-c('Intercept','CCA','Foliose HC', 'Massive zoanthids',
                         'Filamentous algae', 'Digitate OC')
D816_nt<-nice_table(stats_table_D816,   highlight=T, note="* p < .05, ** p < .01, *** p < .001")
D816_nt

Term

β

SE

t

p

95% CI

Intercept

0.764

0.00

145.70

< .001***

[0.75, 0.77]

CCA

-0.409

0.12

-3.46

.001**

[-0.64, -0.18]

Foliose HC

-0.052

0.03

-1.53

.125

[-0.12, 0.01]

Massive zoanthids

0.091

0.06

1.63

.104

[-0.02, 0.20]

Filamentous algae

0.042

0.03

1.26

.207

[-0.02, 0.11]

Digitate OC

0.470

0.33

1.43

.154

[-0.18, 1.12]

Note. * p < .05, ** p < .01, *** p < .001

#D1632
D1632_CI<-confint(D1632_mod)[-nrow(confint(D1632_mod)),]
stats_table <- as.data.frame(D1632_mod_sum$coefficients$cond)
stats_table_D1632<-cbind(
  row.names(stats_table),
  round(stats_table,3), D1632_CI[,-3]
)
colnames(stats_table_D1632) <- c(
  "Term", "B", "SE", "t", "p",
  "CI_lower", "CI_upper")
stats_table_D1632$Term<-c('Intercept','CCA')
D1632_nt<-nice_table(stats_table_D1632,   highlight=T, note="* p < .05, ** p < .01, *** p < .001")
D1632_nt

Term

β

SE

t

p

95% CI

Intercept

2.147

0.01

167.63

< .001***

[2.12, 2.17]

CCA

-0.832

0.33

-2.51

.012*

[-1.48, -0.18]

Note. * p < .05, ** p < .01, *** p < .001

#D3264
D3264_CI<-confint(D3264_mod)[-nrow(confint(D3264_mod)),]
stats_table <- as.data.frame(D3264_mod_sum$coefficients$cond)
stats_table_D3264<-cbind(
  row.names(stats_table),
  round(stats_table,3), D3264_CI[,-3]
)
colnames(stats_table_D3264) <- c(
  "Term", "B", "SE", "t", "p",
  "CI_lower", "CI_upper")
stats_table_D3264$Term<-c('Intercept', 'CCA', 'Massive HC')
D3264_nt<-nice_table(stats_table_D3264,   highlight=T, note="* p < .05, ** p < .01, *** p < .001")
D3264_nt

Term

β

SE

t

p

95% CI

Intercept

2.224

0.02

122.72

< .001***

[2.19, 2.26]

CCA

-1.343

0.38

-3.50

< .001***

[-2.10, -0.59]

Massive HC

0.439

0.16

2.82

.005**

[0.13, 0.74]

Note. * p < .05, ** p < .01, *** p < .001

#Sq
VRSD_CI<-confint(VRSD_mod)[-nrow(confint(VRSD_mod)),]
stats_table <- as.data.frame(VRSD_mod_sum$coefficients$cond)
VRSD_CI<-as.data.frame(VRSD_CI)
stats_table_VRSD<-cbind(
  row.names(stats_table),
  round(stats_table,3), t(VRSD_CI[-3,])
)
colnames(stats_table_VRSD) <- c(
  "Term", "B", "SE", "t", "p",
  "CI_lower", "CI_upper")
stats_table_VRSD$Term<-c('Intercept')
VRSD_nt<-nice_table(stats_table_VRSD,   highlight=T, note="* p < .05, ** p < .01, *** p < .001")
VRSD_nt

Term

β

SE

t

p

95% CI

Intercept

0.343

0.03

11.23

< .001***

[0.28, 0.40]

Note. * p < .05, ** p < .01, *** p < .001

#Save tables

flextable::save_as_docx(V4_nt, path = "Figures/V4_GLMM.docx")
flextable::save_as_docx(PROC32_nt, path = "Figures/PROC32_GLMM.docx")
flextable::save_as_docx(PLC4_nt, path = "Figures/PLC4_GLMM.docx")
flextable::save_as_docx(PLC16_nt, path = "Figures/PLC16_GLMM.docx")
flextable::save_as_docx(PLC32_nt, path = "Figures/PLC32_GLMM.docx")
flextable::save_as_docx(D12_nt, path = "Figures/D12_GLMM.docx")
flextable::save_as_docx(D24_nt, path = "Figures/D24_GLMM.docx")
flextable::save_as_docx(D816_nt, path = "Figures/D816_GLMM.docx")
flextable::save_as_docx(D1632_nt, path = "Figures/D1632_GLMM.docx")
flextable::save_as_docx(D3264_nt, path = "Figures/D3264_GLMM.docx")
flextable::save_as_docx(VRSD_nt, path = "Figures/Sq_GLMM.docx")

9.3. Examining the GLMMs (pseudo-Rsquare and model diagnistic)

############################Check pseudo R2##########################################

library(MuMIn)
pseudo_r2_table<-as.data.frame(matrix( ncol=2,  nrow=11))
row.names(pseudo_r2_table)<-c('VRM - 4 cm', 
                              'PROC - 32 cm',
                              'PLC - 4 cm',
                              'PLC - 16 cm',
                              'PLC - 32 cm',
                              'D(1 - 2 cm)',
                              'D(2 - 4 cm)',
                              'D(8 - 16 cm)',
                              'D(16 - 32 cm)',
                              'D(32 - 64 cm)',
                              'Sq')
pseudo_r2_table[1,]<-round(r.squaredGLMM(V4_mod), 4) #model with best conditional R2
pseudo_r2_table[2,]<-round(r.squaredGLMM(PROC32_mod)[1,], 4)
pseudo_r2_table[3,]<-round(r.squaredGLMM(PLC4_mod),4)
pseudo_r2_table[4,]<-round(r.squaredGLMM(PLC16_mod)[1,], 4) 
pseudo_r2_table[5,]<-round(r.squaredGLMM(PLC32_mod), 4)#The least explained model
pseudo_r2_table[6,]<-round(r.squaredGLMM(D12_mod)[1,],4)
pseudo_r2_table[7,]<-round(r.squaredGLMM(D24_mod)[1,],4)# The second best-explained model
pseudo_r2_table[8,]<-round(r.squaredGLMM(D816_mod)[1,],4)
pseudo_r2_table[9,]<-round(r.squaredGLMM(D1632_mod),4)
pseudo_r2_table[10,]<-round(r.squaredGLMM(D3264_mod),4)
pseudo_r2_table[11,]<-round(r.squaredGLMM(VRSD_mod),4)
colnames(pseudo_r2_table)<-c('Marginal R squared', 'Constrained R squared')
#Save tables
library(gt)
library(rempsyc)
pseudo_r2_table$Indicator<-row.names(pseudo_r2_table)
pseudo_r2_table<-dplyr::select(pseudo_r2_table, Indicator, `Marginal R squared`, `Constrained R squared`)
pseudo_r2_nice<-gt(pseudo_r2_table)

pseudo_r2_nice
Indicator Marginal R squared Constrained R squared
VRM - 4 cm 0.7648 0.7648
PROC - 32 cm 0.6110 0.6626
PLC - 4 cm 0.4387 0.4387
PLC - 16 cm 0.2758 0.2758
PLC - 32 cm 0.2447 0.2447
D(1 - 2 cm) 0.6515 0.6515
D(2 - 4 cm) 0.8058 0.8058
D(8 - 16 cm) 0.4387 0.4387
D(16 - 32 cm) 0.2281 0.2573
D(32 - 64 cm) 0.4723 0.4723
Sq 0.0000 0.0000
gtsave(pseudo_r2_nice, 'Figures/PseudoR2.docx')
######################Model diagniostics################################

library(DHARMa)
library(SuppDists)


V4_resi<-simulateResiduals(V4_mod)

plot(V4_resi, sub='VRM - 4 cm model diagnosis')

PROC32_resi<-simulateResiduals(PROC32_mod)

plot(PROC32_resi, sub='PROC - 32 cm model diagnosis')

PLC4_resi<-simulateResiduals(PLC4_mod)

plot(PLC4_resi, sub='PLC - 4 cm model diagnosis')

PLC16_resi<-simulateResiduals(PLC16_mod)

plot(PLC16_resi, sub='PLC - 16 cm model diagnosis')

PLC32_resi<-simulateResiduals(PLC32_mod)

plot(PLC32_resi, sub='PLC - 32 cm model diagnosis')

D12_resi<-simulateResiduals(D12_mod)

plot(D12_resi, sub='D (1 - 2 cm) model diagnosis')

D24_resi<-simulateResiduals(D24_mod)

plot(D24_resi, sub='D (2 - 4 cm) model diagnosis')

D816_resi<-simulateResiduals(D816_mod)

plot(D816_resi, sub='D (8 - 16 cm) model diagnosis')

D1632_resi<-simulateResiduals(D1632_mod)

plot(D1632_resi, sub='D (16 - 32 cm) model diagnosis')

D3264_resi<-simulateResiduals(D3264_mod)

plot(D3264_resi, sub='D (32 - 64 cm) model diagnosis')

VRSD_resi<-simulateResiduals(VRSD_mod)

plot(VRSD_resi, sub='Sq model diagnosis')

10. 3D model quality

Showing the reprojection error, ground resolution, and control error of each model.

rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment")
########################################################################
library(readxl)
mod_qual<-read_xlsx('ModelQuality.xlsx')
## New names:
## • `` -> `...1`
mod_qual<-as.data.frame(mod_qual)
row.names(mod_qual)<-mod_qual[,1]
mod_qual<-mod_qual[,-1]

mod_qual$Reproj_Error_mm<-mod_qual$Res_mm*mod_qual$Reproj_Error

#Showing resolutions, reprojection error (both in pixel and mm), and planar area of each quadrat.
head(mod_qual)
##      images Res_mm Reproj_Error  Cont_Error     PLA Reproj_Error_mm
## 82.5   1545  0.858        0.652 0.009873100 23.9740        0.559416
## BT     1621  0.736        0.554 0.000184904 23.0720        0.407744
## CK     1443  0.998        0.531 0.000607361 24.5164        0.529938
## CJ     1473  0.811        0.591 0.000196119 22.6973        0.479301
## CCG    1657  0.679        0.624 0.000299013 23.6267        0.423696
## DaBS   1527  0.812        0.548 0.000218082 23.7451        0.444976
#Showing the ranges of parameters
mod_qual_sum<-data.frame()

for(i in 1:ncol(mod_qual)){
  mod_qual_sum[1,i]<-min(mod_qual[,i])+(max(mod_qual[,i])-min(mod_qual[,i]))/2
  mod_qual_sum[2,i]<-(max(mod_qual[,i])-min(mod_qual[,i]))/2
  mod_qual_sum[3,i]<-mean(mod_qual[,i])
}
colnames(mod_qual_sum)<-colnames(mod_qual)

row.names(mod_qual_sum)<-c('Intermediate', 'HalfRange', 'Mean')

#Showing the ranges
round(mod_qual_sum, 3)
##              images Res_mm Reproj_Error Cont_Error    PLA Reproj_Error_mm
## Intermediate 1394.5  0.875        0.604      0.005 23.792           0.584
## HalfRange     440.5  0.284        0.136      0.005  2.316           0.262
## Mean         1504.6  0.774        0.606      0.001 24.154           0.473

11. Displaying the maps

rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment")
library(maptools)
library(rgdal)
library (sp)
library(raster)
library(ggplot2)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(ggspatial)
library(rgbif)
#library(mapr)
library(marmap)
library(leaflet)
library(ggsn)
library(RgoogleMaps)
library(ggmap)
library(mapproj)
#install.packages("leaflet")
library(leaflet)
library(devtools)

#Import geographical data of quadrats
library(readxl)
coordi<-read_xlsx('Photogrammetric_Plot_Coordinates.xlsx')
coordi<-as.data.frame(coordi)
coordi<-coordi[,-c(4,5,6)]
row.names(coordi)<-coordi[,1]

#Separate the data by regions
KT<-coordi[c(6:10),]
LY<-coordi[c(11:15),]
LD<-coordi[c(1:5),]
EC<-coordi[c(16:20), ]
NC<-coordi[c(21:25),]
###############################################
library(ggmap)
library(ggsn)
#Installing the key of google map
register_google(key = "AIzaSyCccvogSM3f_4VjI7k7OmvWzA2lAKDUwQc")
#Get map of Taiwan
map <- get_map(location = 'Taiwan', zoom = 7, map='satellite')
#Plot the sampling sites on map
taiwan<-ggmap(map)+
  scale_x_continuous(limits = c(119.2, 122.5), expand = c(0, 0)) +
  scale_y_continuous(limits = c(21.5, 25.6), expand = c(0, 0))+
  scalebar(x.min = 119.5, x.max=120, y.min=22, y.max=22.5,
              dist = 30, model = 'WGS84',dist_unit = 'km',
           transform = T,st.size = 2.5,st.dist=0.3,box.color = 'white',
          box.fill = c('white', 'yellow'),
          st.color='white',
          st.bottom = F) +
  labs(x='Longitudes', y='Latitudes')+
  geom_rect(aes(xmin = 120.6, xmax = 121, ymin = 21.8, ymax = 22.2), 
            color = '#fc0f00', 
            fill = NA,lwd=0.8)+
  geom_rect(aes(xmin = 121.49, xmax = 121.64, ymin = 21.95, ymax = 22.14), 
            color = '#F28522', 
            fill = NA,lwd=0.8)+
  geom_rect(aes(xmin = 121.44, xmax = 121.53, ymin = 22.62, ymax = 22.69), 
            color = "#fec700", 
            fill = NA,lwd=0.8)+
  geom_rect(aes(xmin = 121.22, xmax = 121.85, ymin = 22.86, ymax = 24.49), 
            color = "#15e0fa", 
            fill = NA,lwd=0.8)+
geom_rect(aes(xmin = 121.756, xmax = 121.965, ymin = 25, ymax = 25.164), 
          color = "#0c59fe", 
          fill = NA, lwd=0.8)+
  annotate('text', x=120.5, y=21.7, label='Kenting',size=5,
           color='white')+
  annotate('text', x=121.5, y=21.9, label='Lanyu',size=5,
           color='white')+
  annotate('text', x=121.5, y=22.5, label='Ludao',size=5,
           color='white')+
  annotate('text', x=122, y=22.75, label='East Coast',size=5,
           color='white')+
  annotate('text', x=122.1, y=25.28, label='North Coast',size=5,
           color='white')
taiwan

ggsave('Figures/Taiwan_map.jpg', width=1500, height=2250, units='px')

#Plot map of Kenting
map1 <- get_map(location = c(lon = 120.769022, lat = 21.953408),
               zoom = 12, map='satellite')
KT_map<-ggmap(map1)+
  scale_x_continuous(limits = c(120.68, 120.88), expand = c(0, 0))+
  scale_y_continuous(limits = c(21.88, 22.05), expand = c(0, 0))+
  scalebar(x.min = 120.7, x.max=120.75, y.min=21.89, y.max=21.92,
           dist = 2, model = 'WGS84',dist_unit = 'km',
           transform = T,st.size = 2.5,st.dist=0.3,box.color = 'white',
           box.fill = c('white', 'yellow'),
           st.color='white',
           st.bottom = F) +
  labs(x='Longitudes', y='Latitudes')+
  geom_point(data=KT, aes(x=Long, y=Lat), color='#fc0f00',
             size=4)+
  annotate('text', x=120.71, y=21.937, label='Dinbaisha',size=3,
           color='white')+
  annotate('text', x=120.755, y=21.923, label='Leidashih',size=3,
           color='white')+
  annotate('text', x=120.76, y=21.937, label='Houbihu',size=3,
           color='white')+
  annotate('text', x=120.78, y=21.96, label='Tiaoshih',size=3,
           color='white')+
  annotate('text', x=120.85, y=21.985, label='Jialeshuei',size=3,
           color='white')
KT_map

ggsave('Figures/Kenting_map.jpg', width=1500, height=1500, units='px')


#Plot map of Lanyu
map2 <- get_map(location = c(lon = 121.551093, lat = 22.047084),
                zoom = 12, map='satellite')
LY_map<-ggmap(map2)+
  scale_x_continuous(limits = c(121.47, 121.625), expand = c(0, 0))+
  scale_y_continuous(limits = c(21.99, 22.12), expand = c(0, 0))+
  scalebar(x.min = 121.45, x.max=121.5, y.min=22, y.max=22.025,
           dist = 1, model = 'WGS84',dist_unit = 'km',
           transform = T,st.size = 2.5,st.dist=0.3,box.color = 'white',
           box.fill = c('white', 'yellow'),
           st.color='white',
           st.bottom = F) +
  labs(x='Longitudes', y='Latitudes')+
  geom_point(data=LY, aes(x=Long, y=Lat), color='#F28522',
             size=4)+
  annotate('text', x=121.55, y=22.01, label='Green Grassland',size=3,
           color='white')+
  annotate('text', x=121.50, y=22.05, label='Yayo',size=3,
           color='white')+
  annotate('text', x=121.55, y=22.088, label='Hen Rock',size=3,
           color='white')+
  annotate('text', x=121.575, y=22.094, label='Two Lions',size=3,
           color='white')+
  annotate('text', x=121.585, y=22.058, label='Donching',size=3,
           color='white')
LY_map

ggsave('Figures/Lanyu_map.jpg', width=1500, height=1500, units='px')


#Map of Ludao
map3 <- get_map(location = c(lon = 121.490517 , lat = 22.659331),
                zoom = 13, map='satellite')

LD_map<-ggmap(map3)+
  scale_x_continuous(limits = c(121.455, 121.52), expand = c(0, 0))+
  scale_y_continuous(limits = c(22.625, 22.688), expand = c(0, 0))+
  scalebar(x.min = 121.46, x.max=121.48, y.min=22.63, y.max=22.64,
           dist = 1, model = 'WGS84',dist_unit = 'km',
           transform = T,st.size = 2.5,st.dist=0.3,box.color = 'white',
           box.fill = c('white', 'yellow'),
           st.color='white',
           st.bottom = F) +
  labs(x='Longitudes', y='Latitudes')+
  geom_point(data=LD, aes(x=Long, y=Lat), color="#fec700",
             size=4)+
  annotate('text', x=121.488, y=22.637, label='Dabaisha',size=3,
           color='white')+
  annotate('text', x=121.48, y=22.642, label='Guiwan',size=3,
           color='white')+
  annotate('text', x=121.469, y=22.653, label='Shihlang',size=3,
           color='white')+
  annotate('text', x=121.48, y=22.68, label='Chaikou',size=3,
           color='white')+
  annotate('text', x=121.493, y=22.683, label='Gongguanbi',size=3,
           color='white')
LD_map

ggsave('Figures/Ludao_map.jpg', width=1500, height=1500, units='px')

##Map of East Coast
map4 <- get_map(location = c(lon = 121.608442, lat = 23.934540),
                zoom = 8, map='satellite')

EC_map<-ggmap(map4)+
  scale_x_continuous(limits = c(121.177, 121.932), expand = c(0, 0))+
  scale_y_continuous(limits = c(22.8, 24.58), expand = c(0, 0))+
  scalebar(x.min = 121.6, x.max=121.8, y.min=22.9, y.max=23.2,
           dist = 10, model = 'WGS84',dist_unit = 'km',
           transform = T,st.size = 2.5,st.dist=0.3,box.color = 'white',
           box.fill = c('white', 'yellow'),
           st.color='white',
           st.bottom = F) +
  labs(x='Longitudes', y='Latitudes')+
  geom_point(data=EC, aes(x=Long, y=Lat), color="#15e0fa",
             size=4)+
  annotate('text', x=121.35, y=22.88, label='Jiamuziwan',size=3,
           color='white')+
  annotate('text', x=121.5, y=23.1, label='Jihui',size=3,
           color='white')+
  annotate('text', x=121.65, y=23.62, label='Hsinshe',size=3,
           color='white')+
  annotate('text', x=121.62, y=23.7, label='Jiqi',size=3,
           color='white')+
  annotate('text', x=121.7, y=24.5, label='Fenniaolin',size=3,
           color='white')
EC_map

ggsave('Figures/EastCoast_map.jpg', width=1000, height=3000, units='px')
#Map of North Coast

map5 <- get_map(location = c(lon = 121.883732  , lat = 25.112880 ),
                zoom = 12, map='satellite')
NC_map<-ggmap(map5)+
  scale_x_continuous(limits = c(121.78, 121.95), expand = c(0, 0))+
  scale_y_continuous(limits = c(25.05, 25.21), expand = c(0, 0))+
  scalebar(x.min = 121.9, x.max=121.945, y.min=25.2, y.max=25.21,
           dist = 2, model = 'WGS84',dist_unit = 'km',
           transform = T,st.size = 2.5,st.dist=0.35,box.color = 'white',
           box.fill = c('white', 'yellow'),
           st.color='white',
           st.bottom = F) +
  labs(x='Longitudes', y='Latitudes')+
  geom_point(data=NC, aes(x=Long, y=Lat), color="#0c59fe",
             size=4)+
  annotate('text', x=121.81, y=25.15, label='Chaojing',size=3,
           color='white')+
  annotate('text', x=121.895, y=25.128, label='82.5K',size=3,
           color='white')+
  annotate('text', x=121.91, y=25.14, label='Bitou',size=3,
           color='white')+
  annotate('text', x=121.92, y=25.11, label='Longdong',size=3,
           color='white')+
  annotate('text', x=121.9, y=25.075, label='Meiyenshan',size=3,
           color='white')
NC_map

ggsave('Figures/NorthCoast_map.jpg', width=1500, height=1500, units='px')
#################################################################